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1  Executive  Summary 

This  is  the  final  technical  report  for  the  ARPA/AFOSR-sponsored  project  “A  Wavelet- 
Unified  Design  Tool  System  for  Signal  Processing,”  contract  F49620-92-C-0054.  This 
project  explored  the  theory  of  a  broad  new  class  of  mathematical  transforms  (rank 
M  wavelets),  developed  a  production-quality  software  package  for  the  design  and  im¬ 
plementation  of  these  transforms,  and  applied  them  to  concrete  problems  in  signal 
processing  and  communications.  Significant  results  of  the  three-year  project  include: 

•  Complete  theoretical  treatment  of  rank  M  wavelets,  including  parametrization, 
design  and  implementation  techniques.  Original  work  on  the  smoothness  (reg¬ 
ularity)  of  wavelet  systems  and  its  application  to  image  compression. 

•  Development  of  a  broad  and  flexible  software  system,  “WaveTool”,  for  design¬ 
ing  and  implementing  wavelet  algorithms.  After  installation  at  a  number  of 
government,  academic,  and  industrial  beta  sites,  this  software  was  successfully 
turned  into  a  commercial  product. 

•  Application  of  rank  M  wavelets  and  design  techniques  to  broadband  communi¬ 
cations,  particularly  in  multicarrier  modulation.  Wavelet  transforms  developed 
under  this  contract  have  been  implemented  in  a  chipset  product  for  high-bitrate 
last-mile  telecommunications. 

•  Application  of  the  wavelet  transforms  and  design  techniques  to  the  compression 
of  sonar,  seismic,  and  multispectral  image  data,  as  well  as  improvements  to 
wavelet-based  still  image  compression. 

•  First  significant  applications  of  emerging  multiwavelet  techniques  to  signal  and 
image  processing. 

The  report  that  follows  details  both  our  vision  for  the  project  and  its  substantive 
achievements,  in  the  hope  that  others  may  take  the  work  even  further. 


2  Project  Goals 

Wavelet  transforms,  discovered  in  the  late  1980’s,  offer  a  variety  of  useful  properties 
and  important  new  applications.  The  current  project  sprung  from  the  development, 
at  Aware  and  elsewhere,  of  a  broad  new  class  of  transforms  called  “rank  M  wave¬ 
lets.”  These  new  methods  encompass  traditional  discrete  block  transforms  such  as  the 
Discrete  Fourier  Transform  (DFT),  Discrete  Cosine  Transform  and  Walsh-Hadamard 
tranforms,  as  well  as  Daubechies  wavelets  and  the  Lapped  Orthogonal  Transform 
(LOT)  and  Modulated  Lapped  Transform  (MLT)  of  subband  coding.  In  addition,  a 
wide  range  of  entirely  new  data  transforms  are  described  by  the  rank  M  framework, 
with  many  potential  applications.  The  goals  of  the  project  were  threefold: 

•  To  gain  a  deeper  theoretical  understanding  of  the  space  of  rank  M  wavelets,  in 
both  the  discrete  and  continuous  domains; 

•  To  build  a  wavelet  software  toolbox  (known  as  WaveTool)  that  enables  the 
engineer  to  prototype  and  test  the  new  wavelet  methods  at  his  desktop; 

•  To  develop  rank  M  wavelet  approaches  to  new  application  problems,  such  as 
seismic  and  sonar  signal  processing,  and  communications. 

We  achieved  significant  results  in  each  of  these  three  areas.  Initial  work  concen¬ 
trated  on  answering  a  number  of  theoretical  questions  regarding  rank  M  wavelets, 
resulting  in  complete  parametrizations,  closed-form  constructions,  and  an  analysis  of 
the  regularity  (smoothness)  of  these  function  systems.  We  completed  development 
of  the  WaveTool  wavelet  software  toolbox,  a  GUI-based  system  for  UNIX  and  PC 
platforms.  An  engineer  can  use  WaveTool  to  rapidly  design  a  variety  of  wavelet  filter 
banks  and  graphically  assemble  them  into  arbitrary  tree  structures.  He  or  she  can 
then  apply  the  resulting  transform  to  data  in  a  variety  of  formats,  and  work  with  the 
results  both  visually  and  numerically.  The  WaveTool  software  package  was  released 
to  a  number  of  government,  academic,  and  industrial  “beta”  sites,  and,  with  further 
modifications,  grew  into  a  commercial  product  offered  by  Aware.  In  addition  to  the 
effort  put  into  theory  and  software,  we  applied  rank  M  wavelet  methods  to  a  number 
of  application  areas,  including  the  compression  of  multispectral,  seismic,  and  sonar 
data.  A  surprise  development  was  the  application  of  rank  M  wavelets  to  multicarrier 
modulation  for  broadband  last-mile  telecommunications.  Wavelet  filter  designs  from 
the  WaveTool  software  proved  to  be  a  critical  enabler  for  this  “wavelet  multitone” 
technology.  Finally,  we  applied  project  resources  to  the  development  of  new  “multi¬ 
wavelet”  methods  for  signal  and  image  processing  (denoising  and  data  compression); 
this  work  represents  the  first  significant  application  of  multiwavelets,  with  promising 
results. 
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In  the  sections  to  follow  we  give  a  background  on  rank  M  wavelets  and  filter 
banks,  followed  by  sections  devoted  to  the  results  of  the  project  in  the  three  different 
directions  of  theory,  software,  and  applications. 


3  Review  of  the  Rank  M  Wavelet  Framework 

In  this  section  we  review  the  general  framework  provided  by  rank  M  wavelets,  es¬ 
tablishing  connections  with  the  theory  of  multirate  filter  banks  and  with  mulitscale 
analysis  and  wavelets.  This  will  help  set  the  stage  for  the  discussion  of  our  theoretical 
and  applied  work  to  follow. 

Rank  M  wavelet  theory  provides  a  unified  methodology  for  the  design  and  analysis 
of  discrete  transforms;  its  instances  share  the  following  important  properties: 

•  The  transforms  satisfy  a  perfect-reconstruction  (invertibility)  property  while 
yielding  a  critically  sampled  data  representation; 

•  It  is  possible  to  impose  orthogonality  on  the  transform  basis  functions  so  that 
the  wavelet  transform  preserves  signal  energy; 

•  The  transform  basis  functions  have  compact  support  and  are  represented  by 
digital  FIR  filters; 

•  The  class  includes  basis  functions  and  associated  filters  that  correspond  to  over¬ 
lapped  windows; 

•  The  class  includes  multirate  (polyphase)  digital  filters; 

•  The  transforms  describe  a  signal  in  time,  frequency,  and  scale; 

•  They  have  “fast”  computational  algorithms. 

This  class  of  data  transforms  includes  and  generalizes  orthogonal  block  transforms 
such  as  the  FFT,  DCT,  and  Hadamard  transforms,  filter  banks  such  as  the  Lapped 
Orthogonal  Transform,  and  the  Daubechies  wavelets. 

3.1  Rank  M  Wavelet  Matrices  and  Filter  Banks 

The  basic  building  block  of  rank  M  wavelet  theory  is  a  wavelet  matrix  or  filter  bank, 
which  partitions  a  signal  into  M  frequency  channels  (M  is  an  integer  greater  than  or 
equal  to  2).  This  partitioning,  which  is  achieved  by  convolution  with  the  matrix  rows 
followed  by  downsampling,  is  shown  in  Figure  1.  We  insist  on  orthogonality  of  the 
transform,  i.e.  that  the  energy  in  the  transform  domain  is  equal  to  the  signal  energy, 
and  that  we  are  able  to  invert  the  transform  using  the  same  matrix.  If  the  M  x  Mg 
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wavelet  matrix  A  has  entries  ar,k,  so  that  each  row  of  the  matrix  is  a  filter  impulse 
response,  then  the  orthogonality  condition  can  be  stated  as 

k 

In  other  words,  the  matrix  of  filter  taps  is  orthogonal  to  itself  when  overlapped  at 
shifts  which  are  multiples  of  M.  A  key  consideration  here  is  that  the  genus  g  of  the 
matrix  may  be  greater  than  1,  that  is,  the  transform  may  have  nontrivial  overlap 
(when  ^  =  1,  a  wavelet  matrix  is  simply  an  M  x  M  orthogonal  transform,  such  as 
the  DFT).  Each  input  “frame”  of  M  samples  contributes  to  more  than  one  output 
frame,  when  g  >  1.  This  extension  from  block  transforms  to  overlapping  transforms  is 
important,  for  it  enables  the  signal  processor  to  optimize  the  filters  in  the  transform  for 
additional  properties  such  as  high  stopband  attenuation  (better  subchannelization), 
or  regularity  (smoothness)  in  a  multiresolution  tree  structure. 


Figure  1:  Ideal  frequency  responses  of  an  M-band  filter  bank  (wavelet  matrix). 

Table  1  indicates  where  various  conventional  digital  signal  processing  transforms 
fit  into  the  general  rank  M  framework,  as  the  wavelet  matrix  parameters  M  and  g 
vary.  Every  box  in  the  diagram  describes  a  family  of  filters,  some  of  which  generalize 
classical  filters,  others  of  which  are  entirely  new.  The  top  row  includes  the  rank  2 
Daubechies  wavelets  of  compact  support,  while  the  first  column  encompasses  classical 
block  transforms.  The  interior  of  the  table,  where  M  >  2  and  5^  >  1,  includes 
early  examples  like  the  Lapped  Orthogonal  Transform  and  the  Modulated  Lapped 
Transform,  but  also  much  uncharted  territory.  One  of  the  goals  of  this  project  has 
been  to  provide  design  tools  for  creating  wavelet  transforms  of  arbitrary  rank  M  and 
genus  and  pushing  them  into  new  applications. 

Many  transform  application  domains  fit  within  the  context  of  the  arbitrary  rank 
wavelet  theory.  For  example,  rank  2  overlap  3  wavelets  have  been  found  to  be  optimal 
for  image  and  video  compression  [130],  while  the  JPEG  image  compression  algorithm 
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Table  1:  Conventional  digital  signal  processing  transforms  within  the  framework  of 
the  arbitrary  rank  wavelet  theory.  W{M,g)  denotes  the  set  of  wavelet  transforms 
of  rank  M  and  genus  (overlap)  g,  while  Z)5'=Daubechies’  wavelet  transform  with  2g 
coefficients. 
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is  based  on  a  rank  8  overlap  1  DOT  block  transform  [75].  The  MPEG  Layer  II  audio 
compression  standard  [5]  employs  a  rank  32  genus  16  wavelet  transform.  In  error 
control  coding  [13],  Reed-Solomon  codes  are  based  on  the  overlap  1  (finite  field)  FFT 
and  Hadamard  block  codes  are  based  on  the  overlap  1  Hadamard  matrices.  One  of  the 
greatest  successes  of  this  project  has  been  the  application  of  wavelet  matrix  methods 
to  multicarrier  modulation  for  broadband  communications.  While  the  ANSI  TlEl.4 
multicarrier  standard  for  Asymmetric  Digital  Subscriber  Line  (ADSL)  transmission 
[6]  specifies  the  use  of  a  rank  512  genus  1  Discrete  Fourier  Transform,  Aware’s  sub¬ 
missions  for  the  next-generation  Very  high  bitrate  Digital  Subscriber  Line  (VDSL) 
standard  are  based  on  a  rank  128  genus  6  wavelet  transform  [109],  [89],  [47]. 

Every  rank  M  wavelet  matrix  can  be  realized  as  an  Af-band  multirate  filter  bank 
[116],  with  each  row  of  the  matrix  corresponding  to  a  filter  in  the  filter  bank.  A 
multirate  filter  bank  works  as  shown  in  Figure  2:  a  signal  is  passed  through  a  bank  of 
M  different  filters  (and  these  filters  usually  partition  the  frequency  spectrum),  and  the 
M  results  are  downsampled  (or  decimated)  by  a  factor  of  M,  so  that  while  each  output 
Vk  (called  a  subband)  has  a  sampling  rate  ^  that  of  the  original  signal,  the  totality 
of  the  outputs  retains  the  critical  sampling  rate.  This  filtering-and-downsampling 
operation  is  the  fundamental  building  block  of  a  rank  M  wavelet  transform.  The 
transform  or  subband  domain  provides  a  new  representation  of  the  signal,  in  which 
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processing  such  as  quantization  (for  compression)  or  thresholding  (for  denoising)  or 
peak-picking  (for  feature  extraction)  may  take  place.  In  order  to  invert  the  transform, 
the  subbands  are  upsampled  (interpolated)  by  a  factor  of  M,  passed  through  the 
(time-reversed)  filters  of  the  wavelet  matrix,  and  summed  to  yield  the  reconstructed 
signal.  If  no  processing  was  done  in  the  transform  domain,  condition  (1)  ensures  that 
the  reconstructed  signal  is  identical  to  the  input,  up  to  a  delay. 


Figure  2:  An  M-band  multirate  filter  bank  associated  with  a  wavelet  matrix 

^Ve  consider  briefly  one  of  the  benefits  of  allowing  arbitrary  overlap  or  genus  in 
a  rank  M  wavelet  transform.  The  frequency  responses  Ak{oj)  of  the  M  filters  in  a 
wavelet  matrix  filter  bank  are  designed  to  uniformly  partition  the  frequency  spectrum 
into  M  pieces,  or  subbands.  If  the  filters  were  ideal  “brickwall”  filters,  this  frequency 
partitioning  would  look  like  Figure  1.  In  practice,  each  wavelet  filter  is  an  FIR,  filter 
and  so  it  cannot  offer  such  ideal  sub  channelization.  Consider  a  widely  used  wavelet 
matrix  of  genus  1,  the  Discrete  Cosine  Transform  (DCT).  The  DCT  is  widely  used 
as  a  real-valued  time-to-frequency  transform,  with  the  rank  of  the  transform  giving 
the  number  of  frequency  bins.  Independent  of  the  rank  M,  the  DCT  s  frequency  bins 
are  non-ideal  to  the  same  extent,  with  sidelobes  that  are  13  dB  below  the  main  lobe 
of  each  component  filter.  This  is  shown  in  Figure  3.  By  moving  to  a  higher  genus 
or  overlapped  rank  M  wavelet  transform,  it  is  possible  to  trade  duration  in  time  for 
better  subchannelization,  while  maintaining  orthogonality  of  the  transform.  Figure  4 
shows  the  frequency  responses  of  such  an  overlapped  wavelet  matrix;  the  example  has 
rank  16  and  genus  4,  and  was  constructed  with  Aware’s  WaveTool  software.  Each 
tone  now  has  sidelobes  that  lie  35  dB  below  the  main  lobe,  providing  a  superior 
approximation  to  the  ideal  filter  bank  of  Figure  1. 
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Figure  3:  Frequency  responses  of  the  filters  in  a  rank  16  Discrete  Cosine  Transform. 


Figure  4:  Frequency  responses  of  the  filters  in  a  rank  16  genus  4  wavelet  matrix 
constructed  using  the  WaveTool  software. 
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As  we  will  see  below,  the  possibility  of  using  overlapped  transforms  also  makes 
possible  the  construction  of  wavelets  (multiresolution  analyses)  with  both  compact 
support  and  underlying  smoothness.  This  was  originally  seen  by  Daubechies  [23],  and 
the  extension  to  arbitrary  values  of  M  was  one  of  the  goals  and  accomplishments  of 
this  research  project. 

3.2  Rank  M  Wavelets  and  Multiscale  Decompositions 

The  essence  of  the  wavelet  idea  is  to  assemble  such  building  blocks  into  a  tree  struc¬ 
ture,  providing  a  decomposition  of  a  signal  into  multiple  resolutions,  with  nonuniform 
frequency  partitioning.  The  original  such  wavelet  decomposition  was  the  Mallat  tree 
[64],  which  iterates  on  the  lowpass  output  of  a  rank  two  wavelet  building  block  to 
create  an  octave-band  decomposition  (Figure  5).  Such  an  octave-band  decomposi¬ 
tion  matches  the  human  visual  system,  and  so  is  well-suited  for  transform  coding  of 
images. 


I - 

_  I - -  HLL 

—  I -  HL 

-  H 

Figure  5:  3-level  Mallat  tree  made  up  of  rank  2  wavelets 

For  this  rank  2  tree,  Mallat  [64]  and  Daubechies  [23]  established  a  precise  mathe¬ 
matical  connection  between  the  operation  of  the  two-channel  wavelet  filter  bank  and 
an  underlying  continuous-time  basis  of  L^(R).  In  particular,  they  found  that  the  im¬ 
position  of  additional  linear  conditions  on  the  filter  bank  or  wavelet  matrix  enabled  the 
Mallat  iteration  to  converge  to  a  basis  function  for  continuous-time.  Not  only  is  this 
condition  necessary  for  the  infinite  iteration  to  make  sense;  it  even  makes  a  difference 
for  a  wavelet  tree  decomposition  of  a  finite  size  [87].  Figure  6  compares  the  overall 
impulse  response  that  results  from  iterating  a  Smith-Barnwell  8-tap  orthogonal  non¬ 
wavelet  filter  eight  times  on  the  lowpass  output  with  the  impulse  response  of  the  same 
iteration  using  a  Daubechies  8-tap  wavelet  filter.  The  iterated  Smith-Barnwell  filter 
displays  a  fractal  behavior  which  makes  it  inappropriate  for  use  in  a  tree-structured 
decomposition.  Further  “regularity”  or  “vanishing  moment”  constraints  on  the  low- 
pass  filters  lead  to  smoother  continuous-time  bases,  and  smoother  finite  iterates  (the 
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Daubechies  8-tap  filter  used  in  Figure  6  has  regularity  order  4).  Whereas  overlapping 
filter  banks  (genus  g  >  1  were  originally  used  to  obtain  better  frequency  localization 
than  block  transforms,  in  the  wavelet  setting  overlap  can  be  used  to  obtain  smooth¬ 
ness  of  the  iterated  filter  bank.  We  have  generalized  the  vanishing  moment  conditions 
to  the  rank  M  case  [100],  [39]  and  developed  explicit  formulae  for  the  construction  of 
rank  M  wavelets  satisfying  A^-th  order  regularity  conditions.  Furthermore  we  have 
discovered  a  set  of  results  relating  these  discrete-time  conditions  to  the  smoothness 
(Sobolev  exponent)  of  the  corresponding  continuous-time  basis  functions.  In  partic¬ 
ular,  we  have  described  the  asymptotic  regularity  of  several  different  constructions 
for  rank  M  wavelet  bases  as  the  filter  length  (support)  increases.  These  issues  are 
discussed  in  detail  in  section  4.3.  Given  a  means  for  measuring  smoothness,  it  can  be 
turned  into  a  criterion  for  wavelet  filter  design,  a  technique  we  have  exploited  in  the 
rank  2  and  rank  M  cases  (sections  4.4.3  -  4.4.6). 


Figure  6:  Top:  Impulse  response  of  upsampling  and  convolution  with  the  Smith- 
Barnwell  8-tap  filter,  iterated  eight  times.  Bottom:  Impulse  response  of  upsampling 
and  convolution  with  the  Daubechies  8-tap  wavelet  filter,  iterated  eight  times. 

The  notion  of  assembling  wavelet  filter  banks  into  arbitrary  tree  structures,  not 
just  Mallat  trees  which  iterate  on  the  lowpass  alone,  is  an  essential  part  of  our  ap¬ 
proach.  Arbitrary  trees  assembled  out  of  rank  M  wavelets  enable  a  wide  range  of 
nonuniform  frequency  decompositions.  Aware’s  CD-quality  real-time  audio  compres¬ 
sion  algorithm  [98]  relies  on  a  two-level  tree  beginning  with  a  rank  8  filter  bank,  and 
iterating  on  5  of  the  8  outputs  with  wavelets  of  varying  rank;  this  tree  is  depicted  in 
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Figure  7.  This  higher  rank  wavelet  decomposition  produces  a  frequency /scale  par¬ 
titioning  which  matches  the  human  ear  -  a  distinctly  different  decomposition  from 
what  one  would  use  for  natural  images,  or  seismic  data.  Each  application  will  have 
its  own  best  decomposition.  Understanding  and  applying  rank  M  multiresolution 
decompositions,  as  well  as  building  a  set  of  software  tools  to  help  propagate  these 
techniques,  have  been  the  goals  for  the  Higher  Rank  Wavelet  project. 


Figure  7:  Tree  used  for  decomposition  of  audio  signals  in  Aware’s  real-time  CD-quality 
compression  algorithm. 


4  Theoretical  Results 

In  this  first  of  three  sections  reviewing  the  technical  results  of  the  Higher  Rank  Wave¬ 
let  Toolbox  project,  we  examine  the  theoretical  results  on  higher  rank  wavelets  that 
were  obtained.  The  structure  and  parametrization  of  rank  M  discrete  wavelets  are 
discussed,  as  well  as  means  for  measuring  the  regularity  of  the  associated  continuous¬ 
time  bases.  We  then  cover  a  number  of  different  techniques  that  have  been  developed 
for  rank  M  wavelet  filter  design,  or  in  mathematical  terms,  the  construction  and 
optimization  of  finite  rank  M  wavelet  sequences  with  desirable  properties  such  as 
subchannelization  or  smoothness  of  the  associated  continuous-time  wavelet  basis. 
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4.1  Structure  of  Rank  M  Wavelet  Systems 

A  discrete  rank  M  wavelet  system  [35,  46,  125, 131]  has  one  scaling  sequence  (lowpass 
filter)  and  M  —  1  wavelet  sequences  (bandpass  and  highpass  filters).  A  real  rank  M 
wavelet  system  is  given  by  an  M  x  K  wavelet  matrix  A  whose  entries  as,k  satisfy 


'Yl^s,kO-s',k+Ml  =  M5s,s'^0,l  ) 
k 

(2) 

and  —  MSs,o  • 

(3) 

k 


The  rows  of  A  are  orthogonal  under  shifts  of  M.  In  contrast  to  the  rank  2  case, 
where  the  single  wavelet  sequence  is  determined  by  the  scaling  sequence,  the  rank 
M  case  has  considerable  freedom  in  the  choice  of  the  M  —  1  wavelet  sequences.  The 
construction  of  rank  M  wavelets  can  be  broken  into  two  steps:  the  design  of  a  lowpass 
filter  (scaling  sequence  ao,k)  and  the  construction  of  a  full  rank  M  wavelet  matrix, 
given  the  scaling  sequence  and  other  information.  We  describe  several  methods  for 
the  design  of  wavelet  lowpass  filters  in  Section  4.4,  and  approaches  to  wavelet  matrix 
construction  in  Section  4.5. 

4.1.1  Wavelet  Bases 

As  observed  in  [46],  [96],  the  wavelet  matrix  describes  a  basis  for  the  space  of  functions 
on  Z.  Specifically,  a  discrete  function  f{k)  may  be  expanded 

M~1  oo 

m)  =  E  E  Cs,ias,Ml+k  (4) 

5=0  — OO 

where  ^ 

=  77  •  (5) 

k 

This  discrete  basis  property  is  equivalent  to  the  orthogonality  condition  (2)  -  the 
wavelet  matrix  provides  a  set  of  overlapping  basis  functions  for  ^^(Z).  The  additional 
“low-pass”  condition  (3)  for  wavelets  confers  the  ability  to  develop  orthonormal  bases 
of  L\R). 

It  is  natural  to  take  these  rank  M  discrete  wavelet  systems  and  seek  to  construct 
compactly  supported  rank  M  wavelet  orthonormal  bases  of  L^(R);  initial  steps  in  this 
direction  have  been  taken  by  Gopinath  and  Burrus  [35].  One  begins  by  considering 
the  rank  M  scaling  function  f  associated  with  the  scaling  sequence  ao,fc,  which  is  the 
solution  to  the  dilation  equation 

m  ao^k(f>{Mx  -  k)  .  (6) 

k 
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Thus  there  is  a  one-one  correspondence  between  the  scaling  sequences  discussed  here 
and  scaling  functions  on  R.  The  dilation  equation  (6)  has  a  unique  compactly  sup¬ 
ported  solution  in  its  Fourier  transform  0  satisfies 

CO 

m  =  TlA,{UM>). 

j~l 

The  corresponding  family  of  wavelets  V’r.j.A:  (1  <  r  <  M  -  1)  is  defined  from  <{>  by 

=  Y^O‘r,k4>iMx  -  k) 

k 

-  k) . 

Lawton  [60]  proved  that  the  collection  form  a  tight  frame  for  L^(R). 

The  explicit  formulae  for  scaling  sequences  developed  below  enable  one  to  con¬ 
struct  rank  M  wavelet  bases  with  arbitrary  smoothness  (as  measured  by  Sobolev 
differentiability),  as  discussed  in  section  4.3. 

4.1.2  The  Polyphase  Matrix  Representation 

Matrices  satisfying  the  orthogonality  condition  (2)  have  been  extensively  studied  in 
the  signal  processing  literature  [113,  116,  118]  under  the  name  “M-band  paraunitary 
perfect  reconstruction  filter  banks.”  Engineers  often  restate  (2)  in  the  2^-transform 
domain,  as  follows:  form  the  M  x  M  polyphase  matrix  H{z)  with  polynomial  entries 

I 


Observe  that  H  and  A  are  related  by 


H(z)  =  Ao  +  zAi  +  . . .  +  . 

H(2:)  is  said  to  be  paraunitary  if  ^^'H.{z)  is  unitary  on  the  unit  circle: 

=  Ml  ,  for  k|  -  1  .  (7) 

Comparison  of  coefficients  of  powers  of  2:  shows  that  the  paraunitarity  of  H  is  equiva¬ 
lent  to  the  orthogonality  under  shifts  (2)  of  the  wavelet  system.  Paraunitary  matrices 
and  polyphase  factorizations  have  been  investigated  in  great  detail  [116].  We  impose 
the  additional  linear  condition  (3)  to  form  a  wavelet  matrix  —  this  amounts  to  the 
requirement  that  the  matrix  Ji(z)\z=i  =  Hq  be  a  Haar  wavelet  matrix.  This  proves 
essential  for  the  development  of  orthonormal  bases  of  T^(R)  (cf.  [23,  35]). 
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The  orthogonality  condition  (2)  can  also  be  stated  in  the  Fourier  domain  -  for 
each  of  the  sequences  consider  its  “symbol”  or  Fourier  transform 

,  A'-l 

Men  ==  J7  E  ^  0<S<M  . 

k=o 


Then  (2)  is  equivalent  to 

M-i  _ 

E  =  ^s,s' 

m=0 

The  paraunitarity  (7)  of  H(2:)  implies  that 

Af-l 

EU.COPsi. 

s=0 


(8) 

(9) 


4.2  Vanishing  Moments  and  Polynomial  Approximation  for 
Rank  M  Wavelets 

The  notion  of  vanishing  moments  of  a  wavelet  systems  was  originally  introduced  by 
Daubechies  in  the  rank  2  case  [23]  as  a  means  of  ensuring  differentiability  of  the 
continuous-time  basis  functions.  Since  then,  vanishing  moments  have  come  to  be 
seen  as  important  for  discrete-time  decompositions  as  well,  and  to  offer  important 
information  on  the  power  with  which  wavelet  representations  can  approximate  poly¬ 
nomials  [101].  As  in  the  rank  2  case,  there  are  several  equivalent  definitions  for  a 
rank  M  wavelet  matrix  A  =  {as,k}  to  have  N  vanishing  moments,  which  we  have 
identified  [100,  39]  in  the  following  theorem: 

Theorem  4.1  A  rank  M  wavelet  system  is  said  to  have  approximation  degree  N  if 
any  of  the  following  equivalent  conditions  hold: 

(i)  [df  d^Yxfr{D)  —  0,  for  n  =  0,  l,...,Af— 1  and  r  =  1,  2,  . . . ,  M  —  1  . 

(ii)  [dfdio)'^Ar{0)  =  0,  for  n  =  0,  1,  . . . ,  A  —  1  and  r  =  1,  2,  . . . ,  M  —  1  . 

(hi)  {d/duYAo  (^)  =0,  /orn  =  0,  1,  . . . ,  A’  -  1;  A;  =  1,  2,  . . . ,  M  -  1  . 

(iv)  -f  MlYao,k+Mi  independent  of  k  for  =  0,  1,  . . . ,  M  —  1  and 

n  =  0,  1,  ...,  A  -  1. 

(v)  Ao{oj)  =  (1  -b  6*“"  4-  . . .  -t-  some  trigonometric 

polynomial  Q(u)). 

Strang  [101]  gives  additional  formulations  of  “degree  A,”  relating  it  to  approxi¬ 
mation  of  order  A  for  functions  on  L^(R),  i.e.  in  the  continous-time  domain. 

These  vanishing  moment  conditions  may  be  used  in  the  design  of  wavelet  lowpass 
filters;  indeed,  the  definition  and  systematic  exploitation  of  rank  M  vanishing  moment 
constraints  is  described  in  Section  4.4  below. 
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4.3  Smoothness  (Regularity)  of  Rank  M  Wavelets 

Once  one  has  constructed  a  compactly  supported  continuous-time  wavelet  basis  from 
a  rank  M  wavelet  matrix,  it  is  natural  to  examine  the  regularity  (differentiability) 
of  the  basis  elements  'tpr{x').  We  have  developed  mathematical  tools  for  measuring 
regularity  in  the  rank  M  setting,  and  explicit  constructions  of  several  families  of 
arbitrarily  differentiable  rank  M  wavelet  bases.  We  describe  our  theoretical  results 
here;  the  constructions  are  described  in  section  4.4  below. 

An  early  theorem  on  wavelets  ([24],  p.l53)  states  that  if  a  wavelet  ^rix)  is  N 
times  continuously  differentiable  and  has  sufficient  decay  to  be  in  L^(R),  then  V'r 
must  vanish  to  order  fV  at  ^  =  0.  This  led  to  the  imposition  of  vanishing  wavelet 
moments  (one  of  the  equivalent  conditions  of  Theorem  4.1)  for  the  construction  of 
wavelet  systems.  Since  the  wavelets  fpr,j,k  s-re  linear  combinations  of  dilates  of  the 
scaling  function  (f>,  the  question  of  their  differentiability  reduces  to  “How  differentiable 
can  we  make  (f>{x)V’  One  standard  way  to  measure  differentiability  of  a  function  in 
L^(R)  is  via  Sobolev  spaces,  which  measure  the  number  of  LP'  derivatives  a  function 

has. 

Given  a  real  number  s,  the  Sobolev  space  Vf  is  defined  by 


Recall  that  for  /  :  R  ^  C,  /  G  =5>  /  €  C“,  for  a  <  s  -  That  is,  if  s  >  n  +  i, 
then  /  has  n  continuous  derivatives.  We  define  the  Sobolev  smoothness  of  a  function 
(/>  G  T^(R)  to  be 

s[(j))  :=  sup  {s  :  ^  G  'H®}  • 

Daubechies  and  others  [23],  [31],  [120],  [122]  developed  theoretical  tools  for  mea¬ 
suring  the  Sobolev  exponent  of  a  rank  2  scaling  function,  and  found  that  the  smooth¬ 
ness  of  a  minimal  length  rank  2  wavelet  system  with  approximation  degree  N  {N 
vanishing  moments)  grows  linearly  with  N .  Specifically,  if  ^jv  denotes  the  scaling 
function  for  such  a  system, 

s(</)iv)  ~  .2075iV  . 

We  set  out  to  generalize  this  result  to  the  rank  M  case,  and  have  found  surprising 
results  [52]. 

As  in  the  rank  2  case,  the  Sobolev  regularity  of  the  scaling  function  (f)  is  deter¬ 
mined  by  the  maximum  eigenvalue  of  a  finite-dimensional  operator  associated  with 
the  scaling  sequence  ao,k-  In  particular,  consider  the  modulus  squared  of  the  lowpass 
frequency  response,  denoted  by  P: 


P{u)  =  \Ao{Lo)f 


1  +  -I- . . .  + 

M 


2N 

R{w)  . 


(10) 
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The  trigonometric  polynomial  R{u!)  is  defined  as  the  remainder  after  division  of  P{u>) 
by  N  powers  of  the  Haar  polynomial  squared;  this  remainder  will  play  an  important 
role  in  both  our  measure  of  smoothness  and  in  later  filter  designs.  Once  given  the 
functions  P{uj)  and  R{u)),  we  define  the  associated  transfer  operators  [21]  that  take 
functions  on  [0,7r)  to  functions  on  [0,7r): 


Tpu{oj)  =  P 


u 

M 


mJ  +  ^  {~iir 


u 


a;  +  27r 

M 


+ 


p  P  27r(M  “  +  27r(M  —  1)' 


M 


M 


and 


I  ^  \  „/a;  +  27r\ 


w  +  27r 

M 


+ 


ppf'^  +  27r(M  -  1)^  ^  /"u;  +  27r(M  -  1)' 


M 


M 


Here  we  have  extended  the  domain  of  definition  of  u  by  its  natural  periodization. 

While  they  are  defined  on  the  infinite-dimensional  space  C[0,7r],  these  trans¬ 
fer  operators  leave  certain  finite-dimensional  subspaces  invariant,  and  their  finite¬ 
dimensional  restrictions  supply  enough  information  to  measure  the  Sobolev  regularity 
of  (j).  We  now  sketch  how  this  is  done.  Let  the  trigonometric  polynomials  P  and  R 
be  given  as  in  (10),  expressed  in  the  form 


L 

Y.  Pi 

j=-L 


iC 


IJUJ 


R(e^ 


')=  Y  Pi^ 


lJ(jJ 


(11) 


Moreover  define  the  2L  -f-  1-dimensional  space  of  trigonometric  polynomials 


Sl:= 


u{u>)  = 


5 


the  analogous  2J  -f-  1-dimensional  subspace  Pj,  and  the  distinguished  subspace 

^L,N  ■=  {/  e  Pl  :  f{io)  =  (1  -  cosw)^5r(cj)} 

of  £i,.  We  have  demonstrated  the  following  theorem,  generalizing  the  rank  2  work  of 
[24],  [31],  [120]: 


Theorem  4.2  Suppose  that  the  wavelet  system  associated  with  the  scaling  sequence 
{ayt}  has  degree  N  (the  first  N  wavelet  moments  vanish).  Then 
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(i)  The  subspace  8l  o/C[0,7r]  is  an  invariant  subspace  for  the  operator  Tp,  i.e. 

Tp  :  Sp  Sp  . 

The  operator  Tp  acting  on  Sp  has  the  matrix  representation  [Tp]j  i^  =  PMj-k,  wUh 

respect  to  the  basis  ,  1, . . . ,  where  pi  =  0  if  I  <  —L  or  I  >  L.  This 

matrix  is  pseudo-circulant,  obtained  from  shifts  of  the  sequence  pi  by  M  for  each  row. 


For  example,  if  L  =  Mg  — 

1  for  some  g  E  Z, 

as  is  the 

case  when 

{afc}  has  length  Mg, 

the  matrix  has  the  fo' 

rm 

0 

0 

0 

0 

0 

0 

0 

0 

PL-M+2 

PL  0 

0 

PL  0 

0 

0 

PL-2M+2 

Pl-m  Pl-m+1 

.  .  . 

0 

Tp  = 

P-L 

PM-2-L  Pm-i-l 

PM-L 

.  .  . 

PL 

. 

0 

0  0 

P-L 

Pl-m 

0 

0 

0  p-p 

... 

Pm-2-l 

0 

0 

0 

0 

0 

0 

0 

0 

(ii)  Similarly,  the  subspace  Sj  ofC[0,7T]  is  an  invariant  subspace  for  the  operator  Tp, 
i.e. 

Tr  :  8j  Sj  . 

We  write  Tr  for  the  restriction  ofTp  to  the  subspace  8j;  it  has  the  matrix 
representation  ^  =  pMj-k,  with  respect  to  the  basis  ,  1, . . . , 

Here  pj  =  0  if  j  <  —J  or  j  >  J .  This  matrix  has  a  pseudocirculant  form  similar  to 
that  ofTp  above. 

(iii)  Tplsj^  includes  among  its  eigenvalues  the  numbers  n  =  0,  1,  . . . ,  2N  —  1. 

(iv)  Suppose  f  G  Hp^N  1  f{^)  =  (1  ~ 

~  (^osufTpgico) .  (12) 
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Hence  J-l,n  is  an  invariant  subspace  for  Tp,  and  moreover  if  f  is  an  eigenvector  for 
Tp  in  J-l,n  with  eigenvalue  then 


TRg{u)  =  M^^Xg{uj) , 


so  that  g  is  an  eigenfunction  for  Tp  with  eigenvalue  M^^X. 
(v)  If  giTp)  is  the  largest  eigenvalue  ofTp\p^jj  and  g{Tp) 
exponent  of  f  is 


s{4>)  =  - 


log  g{Tp) 
21ogM  ■ 


<  I,  then  the  Sobolev 


(13) 


(vi)  If  iji,[Tr)  is  the  largest  eigenvalue  ofTpl^j  j^  and  ^{Tr)  < 


s{(f>)  =  N- 


log/^(rfl) 

21ogM  ■ 


(14) 


Thus  we  have  obtained  an  explicit  formula  for  the  Sobolev  smoothness  of  the 
scaling  function  cj)  in  terms  of  the  maximal  eigenvalue  of  the  associated  multiresolution 
operator  Tp  acting  on  a  distinguished  subspace  of  the  finite-dimensional  space  £n-i, 
and  a  corresponding  formula  in  terms  of  the  maximal  eigenvalue  of  Tr.  These  formulae 
can  be  used  both  to  measure  the  Sobolev  smoothness  of  a  scaling  function,  given  the 
scaling  sequence  that  defines  it,  and  as  a  tool  for  use  in  filter  design.  We  explore  the 
first  of  these  uses  immediately,  and  leave  the  second  to  section  4.4. 

Not  only  is  the  smoothness  determined  by  the  eigenvalue  ^{Tr)  or  g{Tp),  it  has 
been  shown  [31],  [52]  that  one  may  estimate  piTp)  just  from  the  values  of  the  fre¬ 
quency  response  i2(cu): 


Theorem  4.3  We  have 

where 


g{TR)  ^  R{ujc) , 


Mtt 

ujc  =  -rr - for  even  M, 

M-f  W 


and 


Uc  —  n  for  odd  M, 
with  the  approximation  improving  as  N  grows  large. 


This  provides  a  technique  for  investigating  the  asymptotic  differentiability  of  rank 
M  minimal  length  degree  N  wavelet  bases  as  the  degree  N  approaches  infinity.  The 
difference  between  the  even  M  and  odd  M  cases  is  critical;  we  have  proven  [52]: 
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Theorem  4.4  Consider  the  minimal  length  trigonometric  polynomials  yielding  N-th 
order  approximation,  with  remainders  defined  in  (16),  (11)  and  (18)  below.  For 
M  =  2  and  4,  the  asymptotic  value  of  o-®  N  grows  large  is 


Rn 


Mtt 


M 


2N 


M  +  lJ  \  2 

Thus,  if  (pN  Is  the  associated  scaling  function. 


y/N 


s{Pn)  ^ 


,  Mtt 
‘  M+1 


2  log  M 


The  growth  in  smoothness  is  approximately  linear  with  N,  with  a  slope  of  1  - 
.2075  when  M  =  2  and\  -  --  0.0362  when  M  =  4. 

For  M  =  3,  the  asymptotic  behavior  of  Rn{o^c)  is  given  by 


1 


and  thus 


Rn{'^)  ~  ’ 

s{(Pn) 


^/N 

log  N 


4  log  M 

The  growth  in  smoothness  with  increasing  filter  length  is  only  logarithmic  for  M  =  3. 


The  logarithmic  growth  when  M  =  3  is  is  due  to  the  fact  that  we  evaluate  R{oj) 
at  the  critical  point  tt  rather  than  at  We  conjecture  that  the  Sobolev  exponent 

of  pN  grows  as 

log  77 


s{Pn) 


4  log  M 


for  all  odd  M .  In  the  case  of  even  M ,  we  conjecture  that  the  Sobolev  exponent  of  ^]v 
grows  linearly,  with  slope 


SM 


This  approaches  zero  quite  rapidly  as  M  increases.  Numerical  measurements  of  the 
Sobolev  exponents  confirm  these  conjectures  for  sequence  lengths  up  to  100,  and 
dilation  factors  M  <  6.  One  implication  of  this  result  is  that  adding  vanishing 
moments  is  not  a  very  effective  means  of  increasing  the  smoothness  of  rank  M  wavelet 
bases  for  odd  M;  improved  approaches  are  explored  in  sections  4.4.1  and  4.4.6. 
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The  asymptotic  growth  described  in  the  theorem  is  shown  in  Figure  8,  which 
compares  the  true  Sobolev  exponents  of  the  minimal  length  scaling  functions  (as 
computed  in  [52])  with  their  predicted  linear  and  logarithmic  asymptotes.  The  dif¬ 
ference  between  the  linear  growth  of  the  cases  M  =  2,  M  =  4  and  the  case  M  =  3 
is  highlighted  in  Figure  9  which  compares  the  graphs  of  scaling  functions  for  these 
three  dilation  factors,  with  N  —  2  and  N  =  12. 


Figure  8:  Asymptotic  growth  in  Sobolev  smoothness  for  minimal  length  degree  N 
scaling  functions,  A"  =  1,  2,  . . . ,  50.  The  left  graph  shows  the  actual  Sobolev  expo¬ 
nents  and  their  linear  asymptotes  for  even  rank  (M  =  2,  M  =  4,  M  =  6),  while  the 
right  graph  shows  the  Sobolev  exponents  and  their  logarithmic  asymptotes  for  odd 
rank  (M  =  3,  M  =  5). 
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M=2,  N=2,  s=1.0 


M=3,  N=2,  s=0.909 


M=4,  N=2,  s=0.854 


Figure  9:  Scaling  functions  for  M  =  2,3,4  and  N  =  2  (top  row)  and  M  —  2,3,4  and 
=  12  (bottom  row). 


4.4  Methods  for  Wavelet  Filter  Design 

In  this  section  we  describe  a  number  of  methods  that  we  have  developed  and  applied 
for  the  design  of  wavelet  lowpass  filters  (scaling  sequences),  including  the  imposition 
of  vanishing  wavelet  moments  (polynomial  approximation),  smoothness  of  the  wave¬ 
let  scaling  function  (both  maximally  smooth  and  frequency-sampling-based  “higher¬ 
smoothness”  lowpass  filters),  and  Nguyen’s  quadratic-constrained  design  approach. 
Additional  results  appear  in  [84] 


4.4.1  Approximation-based  Criteria  (Vanishing  Moments) 


Our  first  constructions  of  rank  M  wavelet  lowpass  filters  were  a  class  that  generalize 
Daubechies  orthogonal  constructions  [23]  from  the  rank  2  case.  As  in  [23],  these  con¬ 
structions  take  off  from  the  characterization  of  vanishing  moments  given  in  Theorem 
4.1.  Recall  from  that  theorem  that  the  frequency  response  of  an  orthogonal  wavelet 
lowpass  filter  with  approximation  of  degree  N  {N  vanishing  wavelet  moments)  must 
have  the  form 

Aoioj)  =  (1  -b  6*'“'  4- . . .  4-  . 

Note  that  1  -f  4-  •  •  •  4-  is  the  Fourier  transform  of  the  Haar  scaling  sequence 

{1,1,...,!}  of  length  M.  In  order  to  characterize  minimal  length  rank  M  wavelet 
systems  with  N  vanishing  moments,  we  shall  describe  the  autocorrelation  |Ao(u;)|^ 
and  then  perform  a  spectral  factorization  to  obtain  Aq-  As  in  the  previous  section, 
for  a  wavelet  system  with  N  vanishing  moments,  we  write 


P{u)  =  |Ao(a;)p 


1  -b  4- . . .  4- 

M 


2N 

RM 


(15) 


where  we  have  factored  the  square  modulus  P{uj)  into  two  parts,  N  powers  of  the 
modulus  squared  Haar  polynomial  and  a  term  R{oj)  =  |(5(u;)|^  designed  to  restore 
orthogonality  to  the  overall  symbol. 

We  have  obtained  [39]  a  complete  closed-form  expression  for  the  minimal  length 
“remainder”  term  RNi<^)  for  arbitrary  values  of  M  and  N: 


N-l 

Rn{<^)  =  X)  -  cosu;)"  ,  (16) 

n=0 

where 


n.  — 


E 


'  Mi-1 


—  krj 


X 


)(l-cos7r)-^^^ 


(17) 
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(18) 


for  M  even,  with  Mi  =  y-  ^  odd, 


f  Ml 

r»=  E  n 

fcl+fc2+—+fcMi  ="  1 


2iV  +  -  1 

2Ar-  1 


1  —  cos 


2'!Tm\ 


— 


M 


where  Mi  =  This  completely  describes  the  deterministic  autocorrelation  Piv(w) 

of  the  minimal  length  rank  M  wavelet  lowpass  filter  with  regularity  order  N .  For  low 
values  of  M,  the  spectral  factor  Ao(tu)  may  be  explicitly  computed.  For  example,  the 
M  =  4,  N  =  2  generalization  of  Daubechies  4-coefficient  scaling  sequence  is 


The  general  (non-minimal  length)  autocorrelation  with  degree  N  approximation 
is  of  the  form 

P{u)  =  Pn{^)  +  P(^)  (20) 


where 


and 


P/v(ai) 


1  q.  -p  . . .  4- 

M 


2N 

Rn{u>) 


P{<jj)  =  TnCosnu  , 

n^Ml 


all  subject  to 


P{ijj)  >0  for  w  G  [0,  tt]  . 


Any  rank  M  wavelet  lowpass  filter  with  N  vanishing  moments  can  be  obtained  as  a 
spectral  factor  of  such  a  P(w),  and  any  such  P{lo)  will  have  finitely  many  spectral 
factors.  These  general  parametrizations  [39]  (see  also  [24],  [126])  can  be  used  to 
describe  various  families  of  orthogonal  wavelets,  including  the  smoother  families  of 
section  4.4.6,  the  “coiflets”  of  [24],  and  the  smoother  wavelets  of  [123]. 


4.4.2  Lagrange  M-th  band  filters  and  M-band  wavelets 

We  have  also  found  an  independent  explicit  formula  for  the  minimal  length  modulus- 
squared  frequency  response  Pn{(jj),  based  on  M-th  band  filters  and  Lagrange  inter¬ 
polation.  This  approach  was  introduced  in  the  2-band  case  [10,  101]  to  connect  the 
maximally  flat  filters  of  Herrmann  [53]  with  their  spectral  factors,  the  Daubechies 
wavelets  [23].  Herrmann’s  construction  can  be  generalized  to  the  M-th  band  case, 
yielding  the  minimal  length  M-th  band  lowpass  filter  with  2Af-th  order  flatness  at 
a;  =  0.  The  general  arbitrary-length  solution  follows  by  the  methods  of  the  previous 


24 


section.  Any  lowpass  filter  satisfying  the  orthogonality  condtion  (2)  or  (8)  will  be  a 
spectral  factor  of  an  M-th  band  filter;  in  fact  (8)  simply  states  that  |Ao(u;)p  is  an 
M-th  band  filter.  In  particular,  rank  M  wavelet  lowpass  filters  with  N  vanishing 
moments  will  be  spectral  factors  of  the  M-th  band  filter  H{oj)  we  find  below. 

Recall  that  an  M-th  band  filter  [67,  115]  is  a  symmetric  FIR  filter  whose  impulse 
response  h[n]  satisfies 

'•[^"1 = { J:  otherwL. 

When  used  in  an  M-fold  interpolator,  these  filters  have  the  useful  property  of  preserv¬ 
ing  the  existing  signal  samples,  while  inserting  new  values  in  between.  Conversely, 
M-fold  decimation  of  h  results  in  the  unit  impulse  response.  When  designing  FIR 
M-th  band  filters,  it  is  desirable  to  approximate  the  ideal  lowpass  filter.  When  the 
M-th  band  filter  is  used  in  the  construction  of  an  M-band  wavelet  lowpass  filter,  it  is 
also  desirable  to  have  some  regularity  -  flatness  of  the  frequency  response  at  a;  =  0. 
These  are  not  conflicting  requirements,  for  flatness  at  w  =  0  together  with  i7(0)  =  1 
forces  lowpass  characteristics  on 

The  M-th  band  condition  (21)  can  be  restated  in  the  frequency  domain: 

+  +  +  (22) 

We  find  a  minimal-length  lowpass  filter  which  satisfies  this  condition,  and  van¬ 

ishes  to  order  2N  at  lo  =  0: 

+  .  (23) 

We  do  so  by  requiring  flatness  of  order  2N  at  the  “test  frequencies”  2TTmfM\ 

H{u)  -f  27rm/M)  =  0{\u>'^^)  ,  1  <  m  <  M  —  1  .  (24) 

One  solves  for  these  conditions  by  constructing  each  of  the  polyphase  components  of 
Hiuj)  as  a  Lagrange  interpolate  to  on  the  coset 

{xn  =  Mn  +  k,  n  =  —  A,  1  —  A,  . . . ,  A  —  1}  . 

Our  result,  proven  in  detail  in  [38],  is 

Theorem  4.5  The  minimal  length  M-th  band  lowpass  filter  with  2N-th  order  flatness 
is  given  by 

1  o  M-lN-l  N-l  TufJxU 

=  S  +  5?  S  S  n  ^  cos(Mn  +  k)^  .  (25) 

k=\  n=0  l=-N  ‘ 

In  fact,  H{(jj)  is  just  the  filter  P/v(w)  of  the  previous  section. 
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Example  {N  —  2)  : 


1  M— 1 

=  (26) 

The  stopband  attenuation  of  these  filters  increases  with  N  (for  fixed  M),  as  can  be 
seen  in  the  4-th  band  examples  shown  in  Figure  10.  The  M-th  band  filters  given  by 
formula  (25)  will  not  approximate  the  ideal  lowpass  filter  as  well  as  those  designed  in 
[67]  using  numerical  optimization.  However,  the  filters  given  here  offer  the  advantages 
of  explicit  formulae  and  rational  coefficients,  as  well  as  the  flatness  at  =  0  which  is 
essential  for  the  construction  of  wavelets. 


Figure  10:  Frequency  responses  in  dB  of  minimal  length  4-th  band  maxflat  filters, 
iV  =  1,  2,  4,  8. 

4.4.3  Maximal  Smoothness  Filter  Design  Criteria 

The  results  of  section  4.3  lead  to  a  new  filter  design  criterion  -  we  can  design  FIR 
lowpass  filters  which  lead  to  maximally  smooth  scaling  functions  (infinitely  iterated 
lowpass  filters),  as  opposed  to  traditional  criteria  such  as  f  or  £°°  (Chebyshev)  error. 
Since  the  smoothness  of  finite  iterates  of  the  lowpass  filter  qualitatively  corresponds 
to  the  smoothness  of  the  scaling  function  [87],  we  can  use  the  mathematical  tools  de¬ 
veloped  for  the  continuous-time  setting  to  help  us  design  wavelet  filters  with  smooth 
finite  iterates.  This  is  important  for  applications  of  wavelets  to  image  coding  because 
the  artifacts  resulting  from  lossy  quantization  reflect  the  shape  of  the  underlying  basis 
function.  If  we  can  optimize  the  smoothness  of  the  discrete  wavelet  basis  functions  re¬ 
sulting  from  iteration  of  the  lowpass  filter,  then  we  should  achieve  superior  perceptual 
quality  for  images  compressed  using  these  maximally  smooth  wavelets. 

The  algorithm  for  measuring  the  smoothness  of  a  given  lowpass  filter  is  the  fol¬ 
lowing; 
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•  Form  the  deterministic  autocorrelation  sequence  p[j]  (time-domain  form  of  F(u;)) 
and  the  matrix  Tp  of  Theorem  4.2. 

•  Find  the  maximal  eigenvalue  p{Tp)  of  this  matrix,  other  than  the  known  eigen¬ 
values  n  =  0,  1,  . . . ,  2N  —  1. 

•  The  Sobolev  smoothness  (number  of  derivatives)  of  the  scaling  function 
associated  to  the  lowpass  filter  is 


s{4>)  = 


log  lj.{Tp) 
21ogM  ■ 


An  alternative  algorithm  may  be  employed,  based  on  the  remainder  R{oj)  and  the  re¬ 
duced  transfer  operator  Tp.  Given  either  method  for  measuring  the  wavelet  smooth¬ 
ness  associated  with  a  lowpass  filter,  several  design  approaches  may  be  employed. 
First,  the  smoothness  measure  s  may  be  used  as  the  objective  function  in  an  opti¬ 
mization  to  yield  filters  that  are  maximally  smooth  for  a  given  sequence  length;  this 
technique  is  used  in  sections  4.4.4  and  4.4.1  below.  Secondly,  an  understanding  of 
the  mathematics  of  Sobolev  regularity  may  be  employed  to  intelligently  impose  addi¬ 
tional  zeros  on  the  frequency  response  of  the  wavelet  lowpass  filter,  in  such  a  way  as 
to  increase  the  smoothness  of  the  associated  scaling  function. 


4.4.4  Sobolev-optimal  Quadrature  Mirror  Filters 

We  first  applied  our  Sobolev-optimality  criterion  to  the  design  of  symmetric  near- 
orthogonal  Quadrature  Mirror  Filters  (QMF’s)  [48].  Recall  that  QMF’s  are  two- 
channel  filter  banks  (rank  2  wavelets)  in  which  the  highpass  filter  is  obtained  by  an 
“alternating  flip”  of  the  lowpass  filter.  If  the  lowpass  filter  of  the  pair  is  ao[n],  a 
sequence  of  length  N,  then  the  highpass  filter  is  given  by 

ai[n]  =  (— l)"ao[A’  —  1  —  n]  .  (27) 

This  can  be  represented  in  the  frequency  domain  as 

Ai{oj)  =  e*“'ylo(7r  -  w)  ; 

Ax  is  a  “mirror”  of  Aq.  One  often  combines  this  construction  with  orthogonality  of 
the  filter  bank, 

ao[n]ao[n  -f  2/j  =  do,/  (28) 

n 

the  rank  2  form  of  equation  (2).  It  was  shown  rather  early  in  the  theory  of  multirate 
filter  banks  [94],  [116]  that  the  only  symmetric  two-channel  orthogonal  QMF  pair  is 
the  Haar  filter,  for  which  uq  =  ||,  ||  .  This  filter  is  too  crude  for  most  image  cod¬ 
ing  applications.  However,  symmetry  (linear-phase)  for  filter  banks  is  an  important 
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property  in  wavelet  and  subband  image  coders;  it  enables  the  use  of  nonexpansive 
symmetric  extension  methods  at  signal  boundaries  [95].  One  solution  to  this  problem 
is  to  give  up  the  QMF  condition  and  seek  filter  banks  which  are  perfect-reconstruction 
(biorthogonal),  rather  than  paraunitary  (orthogonal)  [11],  [73],  [117].  This  approach 
led  to  the  well-known  (7,9)-tap  biorthogonal  wavelet  pair  used  in  the  FBI’s  fingerprint 
compression  standard. 

A  different  solution  to  the  problem  of  symmetry  and  orthogonality  was  developed 
by  Adelson  and  Simoncelli  [8];  they  mildly  relaxed  the  orthogonality  condition  (28)  and 
by  doing  so  were  able  to  obtain  symmetric  filters!  The  odd-length  filters  they  designed 
satisfied  the  orthogonality  conditions  (28)  within  l.Oe-3,  and  led  to  near-perfect- 
reconstruction  of  wavelet-transformed  images.  The  use  of  odd-length  filters  is  also 
noteworthy,  for  the  whole-sample  symmetry  of  an  odd-length  filter  preserves  centers 
of  mass,  leading  to  better  retention  of  edge  details  in  a  lossy  subband  compression 
scheme.  Even-length  filters  will  smear  single-pixel  detail  across  two  filter  outputs, 

leading  to  smoothed,  blurry  artifacts. 

Several  constraints  are  imposed  on  the  filters:  we  constrain  the  length. 

E“n  =  l.  P9) 

n 

and  we  impose  one  (and  therefore  two)  vanishing  wavelet  moments: 

n 

=  0  .  (31) 

n 

Observe  that  if  the  filters  were  truly  orthogonal  (i.e.  fully  satisfied  the  conditions 
(28))  then  one  of  the  equations  (30)-(31)  would  be  superfluous.  However,  in  the 
near-orthogonal  case  of  symmetric  QMF’s,  both  equations  are  necessary.  This  leaves 
two  free  parameters  for  a  nine-tap  QMF  (one  parameter  for  the  seven-tap  case).  We 
optimized  to  find  the  filter  with  the  greatest  Sobolev  smoothness  while  satisfying  the 
orthogonality  constraints  (28)  to  within  a  reasonable  epsilon  (less  than  le-03).  This 
amounts  to  the  nonlinearly  constrained  optimization  problem  of  finding  the  solution 
which  minimizes  the  deviation  from  orthogonality  while  maximizing  the  Sobolev  ex- 
ponent. 

This  approach  led  to  a  new  family  of  QMFs;  the  7- and  9-tap  examples  are  given  m 
Table  2  below.  The  9-tap  Sobolev-optimal  QMF  is  quite  distinct  from  the  9-tap  QMF 
of  Adelson-Simoncelli[8];  its  frequency  response  is  plotted  in  Figure  11,  while  pictures 
comparing  the  corresponding  scaling  and  wavelet  functions  are  shown  in  Figure  12. 
Observe  that  our  Sobolev-optimal  design  yielded  a  flatter  frequency  response  than 
the  previous  frequency-sampling  approach  [8].  The  scaling  function  associated  with 
the  Sobolev-optimal  QMF  has  Sobolev  exponent  s  =  1.795,  and  is  visibly  smoother 
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than  the  scaling  function  associated  with  the  Adelson-Simoncelli  QMF,  which  has 
s  =  1.567.  Both  of  these  scaling  functions  are  smoother  than  the  analysis  scaling 
function  of  the  Daubechies  (7,9)-tap  pair,  which  has  s  =  1.410.  The  Sobolev-optimal 
9-tap  QMF  has  a  maximum  deviation  from  orthogonality  of  6.6e  -  04,  superior  to 
the  7.9e  —  04  max  deviation  for  the  Adelson-Simoncelli  9-tap  QMF.  The  new  7-tap 
Sobolev-optimal  QMF  proved  to.be  quite  similar  to  Adelson  and  Simoncelli’s  7-tap 
design. 


n 

ao[n] 

n 

ao[n] 

0 

-7.8125000e-03 

0 

1.9531250e-02 

1 

-7.3069675e-02 

1 

-4.6875000e-02 

2 

3.6136589e-01 

2 

-7.3234012e-02 

3 

8.5324613e-01 

3 

4.0042839e-01 

4 

8.1451231e-01 

Table  2:  7  and  9-tap  Sobolev-optimal  QMF  lowpass  filters. 


Adelson-Simoncelli  9-tap  QMF 


Optimally  smooth  9-tap  QMF 


Figure  11;  Frequency  responses  of  Adelson-Simoncelli  and  optimally  smooth  lowpass 
filters. 

We  applied  the  new  QMF  banks  to  lossy  wavelet-based  image  compression  of 
natural  and  fingerprint  images,  and  found  the  9-tap  QMF  to  be  competitive  with 
the  Adelson-Simoncelli  QMF  in  both  numerical  and  qualitative  senses;  both  of  these 
wavelet  filters  were  superior  to  the  Daubechies  (7,9)-tap  pair  at  compressing  a  finger¬ 
print  image  with  the  FBFs  WSQ  tree.  Further  details  are  given  in  [48]. 
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Scaling  function  for  Adelson-Simoncelii  9-tap  QMF,  s  =  1 .567 


Scaling  function  for  optimally  smooth  9-tap  QMF,  s  =  1.795 


Figure  12:  Comparison  of  Adelson-Simoncelii  and  optimally  smooth  scaling  functions. 


4.4.5  Maximally  Holder-smooth  Wavelet  Filters 

The  maximal-smoothness  approach  was  also  applied  to  the  design  of  2-channel 
orthogonal  wavelets,  but  this  time  using  the  time-domain  notion  of  Holder  regularity, 
as  opposed  to  the  L^-frequency  domain  notion  of  Sobolev  regularity.  We  now  briefly 
summarize  our  approach  and  the  results;  further  details  appear  in  [58]. 

Our  method  utilized  the  autocorrelation  sequence  domain  where  1)  the  dimension 
of  the  parameter  space  can  be  considerably  reduced  and  2)  an  unconstrained  opti¬ 
mization  problem  results  (as  opposed  to  a  constrained  problem  if  one  were  to  work  in 
the  impulse  response  domain).  We  maximized  the  Holder  smoothness  over  the  class 
of  wavelet  autocorrelation  sequences  of  a  given  length. 

Holder  smoothness  is  a  direct  generalization  of  the  notion  of  differentiability  from 
the  natural  to  the  real  numbers.  A  function  /  has  Holder  smoothness  of  order  a  G 
[0,  1)  if,  for  any  t,  /i, 

|/((  +  A)  -  /(t)l  <  c\h\°,  (32) 

where  c  <  oo  is  a  constant  independent  of  t  and  h.  For  a  higher  order  of  Holder 
smoothness  S}j{f)  =  L  +  a,  L  G  N,  a  G  [0,  1),  the  definition  above  is  applied  to 
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the  L-th.  derivative  of  /.  It  can  be  shown  that  a  function  /  possesses  L  continuous 
derivatives  if  and  only  if  it  has  Holder  smoothness  of  order  -Siy(/)  >  L.  Rioul  has 
given  an  algorithm  for  determining  the  Holder  smoothness  [86]  that  can  be  easily 
implemented  and  yields  a  very  accurate  estimate  for  sh.  Given  the  definition  of 
Holder  smoothness,  the  optimization  problem  to  be  solved  can  be  stated  as  follows: 

Maximize  sh{(p)  over  the  coefficients  p  =  \pK,  Pk-2,  •  •  • ,  Pi]'^  of  P{u>) 
according  to  equations  (2)-(3),  and  subject  to 

P{oj)  >  0.  (33) 

In  order  to  compute  the  scaling  filter  (and  not  its  autocorrelation  sequence),  one  must 
find  the  spectral  factor  [24]  of  the  solution  p. 

The  solution  employs  a  general  purpose  optimization  algorithm  requiring  only  a 
function  that  computes  the  order  of  smoothness  sh  for  any  valid  set  of  parameters. 
The  approach  outlined  here  works  as  well  for  any  other  method  of  measuring  smooth¬ 
ness,  such  as  Sobolev  differentiability  ([31],  [48],  [120]);  and  further  explorations  are 
described  in  [59]. 

It  is  obviously  desirable  to  reduce  the  number  of  free  parameters  to  reduce  the 
complexity  of  the  problem.  One  can  do  so  by  means  of  the  previously  quoted  theorem 
relating  wavelet  differentiability  to  vanishing  moments  [24].  For  a  wavelet  scaling 
function  to  have  N  continuous  derivatives,  the  ^-transform  P{z)  must  have  2{N  -j-  1) 
zeros  at  2;  =  —1.  Thus  the  dimension  of  the  parameter  space  can  be  reduced  to 
—  N.  Although  we  do  not  know  the  achievable  order  of  smoothness  for  a  given 
filter  length  /T  -f  1,  we  have  a  lower  bound  by  measuring  the  Holder  smoothness  of 
the  corresponding  Daubechies  filter.  For  example,  the  D6  filter  has  3  free  parameters 
but  Sff  1.08,  so  that  there  is  just  3  —  —  1  =  1  parameter  left  for  optimization. 

The  remaining  free  parameters  are  further  constrained  by  the  nonnegativity  con¬ 
dition  (33).  It  effectively  constrains  the  possible  root  locations  of  P{z)-  Single  roots 
on  the  unit  circle  are  not  admissible;  complex  roots  off  the  unit  circle  have  to  occur 
in  quadruples,  real  roots  must  occur  in  mirror  pairs,  and  roots  on  the  unit  circle  must 
come  in  complex  conjugate  pairs  with  multiplicity  two,  thus  fixing  two,  one,  and  two 
parameters,  respectively.  It  is  important  to  note  that  each  of  the  three  possible  cases 
can  be  expressed  by  equations  that  are  linear  in  the  autocorrelation  coefficients. 

However,  it  is  not  clear  a  priori  which  of  the  different  combinations  of  roots  (e.g., 
in  the  case  of  two  free  parameters  either  one  complex  quadruple  or  a  pair  of  double 
roots  on  the  unit  circle  can  be  chosen)  leads  to  the  maximum  smoothness.  Thus  one 
has  to  try  all  possible  constellations.  Our  algorithm  may  be  summarized  as  follows: 

•  Choose  the  length  K  +  l  of  the  scaling  filter. 

•  Determine  the  number  of  degrees  of  freedom,  i.e.,  find  how  many  vanishing 
moments  can  be  relaxed. 
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Figure  13:  Holder  smoothness  for  Daubechies  (x)  and  optimized  (o)  filters  as  a 
function  of  filter  length. 

•  Determine  all  possible  constellations  that  could  give  rise  to  a  nonnegative  power 
spectrum. 

•  For  each  possible  constellation,  numerically  maximize  the  desired  smoothness 
over  the  parameter  space,  i.e.,  the  ^  -  iV  radii  and  angles.  A  general  purpose 
optimization  routine  is  used.  It  only  requires  a  function  that  computes  the  order 
of  smoothness,  given  the  parameter  values. 

We  carried  out  the  above  method  for  filter  lengths  ranging  from  4  to  30.  When  the 
filter  length  was  10  or  larger,  we  obtained  filters  with  greater  Holder  smoothness  than 
the  Daubechies  orthogonal  wavelets.  For  lengths  4,  6,  and  8,  the  optimization  proce¬ 
dure  did  not  yield  higher  smoothness  than  the  Daubechies’  filters.  Figure  13  shows  the 
the  Holder  exponent  sh  as  a  function  of  /F-f  1  for  the  scaling  functions  corresponding 
to  Daubechies  filters  and  to  the  new  filters.  Beyond  length  8  the  new  method  clearly 
improves  the  smoothness,  with  increasing  success  for  larger  filter  lengths.  A  linear 
least-squares  fit  to  each  of  the  smoothness  plots  found  the  Daubechies  filters  with 
order  less  than  30  to  have  a  Holder  exponent  growing  as  .1418(A" -hi),  while  the  max¬ 
imally  smooth  filters  have  a  Holder  exponent  growth  of  .2024(7^ -hi)-  We  know  [122] 
that  asymptotically  the  Daubechies  filters’  Hblder  exponents  grow  as  .1037(/F  -h  1). 
The  asymptotic  behavior  of  the  maximal  possible  Holder  exponent  for  a  given  filter 
length,  as  the  length  increases,  remains  an  open  question. 

In  Figure  14  we  have  plotted  the  number  of  zeros  at  tt  for  the  Daubechies  and 
maximally  smooth  scaling  filters;  this  is  equal  to  the  number  of  vanishing  moments 
of  the  corresponding  wavelet  system.  Notice  that  while  Daubechies  filters  of  length 
K  -I-  1  have  a  zero  of  order  ^  at  tt,  the  maximally  smooth  filters  have  a  zero  of 
roughly  half  that  order.  In  other  words,  we  have  given  up  half  the  vanishing  moments 
of  the  Daubechies  filter  to  obtain  free  parameters  for  optimizing  the  Holder  regularity. 
Notice  that  this  behavior  is  not  monotonic,  probably  because  several  different  root 
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Figure  14:  Order  of  zero  at  pi  (number  of  vanishing  moments)  for  Daubechies  (x) 
and  maximally  smooth  (o)  filters  as  a  function  of  filter  length. 
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Figure  15:  Scaling  functions  (top)  and  their  4-th  derivatives  (bottom)  corresponding 
to  26-coefficient  Daubechies  and  Holder-optimized  filters. 


constellations  lead  to  near-optimal  filters. 

Consider  the  case  of  a  length  26  {K  =  25)  scaling  filter  that  is  optimized  for  sr- 
The  Daubechies  filter  corresponds  to  smoothness  sr  =  4.005.  After  optimization 
of  three  distinct  complex  quadruples  (off  the  unit  circle)  one  gets  a  scaling  function 
with  smoothness  sr  =  5.06,  one  more  continuous  derivative.  Figure  15  depicts  the 
scaling  functions  and  their  4-th  derivatives  for  each  of  these  filters.  While  the  scaling 
functions  have  a  similar  appearance,  the  higher  degree  of  smoothness  of  the  derivative 
is  quite  apparent.  The  filter  frequency  responses  are  shown  in  Figure  16. 
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Figure  16;  Frequency  responses  in  dB  of  the  26-coef-ficient  Daubechies  (left)  and 
Holder-optimized  (right)  filters. 

4.4.6  A  Frequency-sampling  Approach  to  Smooth  Wavelet  Filters 

In  addition  to  simply  optimizing  for  maximally  smooth  filters,  one  may  use  the  un¬ 
derlying  mathematics  of  wavelet  regularity  to  inform  the  design  of  wavelet  lowpass 
filters,  while  using  the  traditional  technique  of  frequency  sampling.  Frequency  sam¬ 
pling  simply  means  specifying  the  value  of  a  filter  response  at  particular  points  in  the 
frequency  domain;  in  this  case  the  value  specified  will  be  zero. 

Recall  that  the  asymptotic  regularity  estimates  of  section  4.3  were  expressed  in 
terms  of  the  values  of  the  remainder  trigonometric  polynomial  R{ui)  at  the  periodic 
point  tOc  of  the  ergodic  map  r  :  tu  -)•  Mu  mod  27r  and  its  preimages.  Notice  as  well 
that  the  Daubechies  and  generalized  Daubechies  scaling  filters  are  specified  by  zeros 
at  the  preperiodic  points  u  =  2'nklM^  1  <  <  M  —  1  (when  M  =  2,  this  reduces 

to  an  N-th.  order  zero  at  tt).  Given  these  examples,  we  have  been  led  to  define  a 
new  class  of  scaling  functions  and  corresponding  wavelet  systems  that  includes  the 
maximal  vanishing  moments  family  as  a  special  case.  The  new  class  consists  of  those 
wavelet  scaling  filters  (i.e.  sequences  satisfying  (2)  and  (3))  whose  Fourier  transform 
has  specified  zeros  of  specified  orders  at  preperiodic  points  of  the  map  r.  The  lowpass 
filters  are  specified  by  a  frequency  sampling  condition.  We  call  these  PPZ  wavelet 
matrices  (for  PrePeriodic  Zero)  and  they  are  determined  by  the  choice  of  preperiodic 
zeros  at  specified  locations  Ui  of  order  W;  this  collection  of  constraints  is  designated 
generically  by  (wi;  Ni,  u^,  N2',  ul,  Nl).  Note  that  the  sum  condition  (3)  simply 
amounts  to  requiring  a  zero  of  order  1  at  27rA:/M,  1  <  /c  <  M  —  1  for  each  PPZ 
wavelet  matrix.  We  have  investigated  several  PPZ  wavelet  matrix  families  for  the 
cases  M  =  2,  3,  4  and  will  show  how  they  can  be  used  to  improve  the  Sobolev 
regularity  of  a  filter  with  a  given  length. 

For  M  =  2,  the  dyadic  wavelet  case,  we  consider  three  PPZ  wavelet  matrix  families 
with  filter  length  2N  and  zeros  of  the  specified  locations  and  orders: 

/jv  :  (tt,  N)  (Daubechies), 

IIn:  (tt,  A-2;  27r/3,  1;  47t/3,  1) 

IIIn  :  (tt,  A  -  4;  27r/3,  1;  47r/3,  1;  Stt/G,  1;  7^/6,  1) 

The  definitions  of  families  IIn  ^•iid  IIIn  were  motivated  by  the  role  of  R{^)  and 
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R{y)  ill  th®  asymptotic  estimates  of  section  4.3.  These  two  families  can  be  con¬ 
structed  explicitly;  instead  of  using  the  minimal  length  remainder  for  N  van¬ 

ishing  moments,  use  the  more  general  expression 

R{u))  =  i?7v(w)  -I-  (1  —  COS(jj)^  ^  fn  cos  UU)  ,  (34) 

n^Mk 

and  choose  the  free  parameters  to  assert  additional  zeros  at  the  points  cu,-.  For 
example,  the  representative  for  family  IIm  with  filter  length  2N  will  be  determined 
by 

R{io)  =  Rm-i{i^)  -|-  (1  —  coscj)^-^  (fi  costu  -|-  r3Cos3u;)  , 
where  fi  and  are  determined  by  the  conditions 


R 


=  0 


=  0  . 


R{uj)  must  have  a  double  zero  at  ujc  =  ^  because  it  is  associated  with  the  magnitude 
squared  of  the  Fourier  transform.  The  coefficients  of  Rn{<^)  are  predetermined  by 
formulae  (16)-(18),  and  can  be  explicitly  solved  for. 

Both  families  II n  and  1 1  In  have  Sobolev  regularity  growing  faster  with  N  than 
the  .2075A^  growth  of  the  Daubechies  orthogonal  wavelets,  and  the  Sobolev  regularity 
of  the  family  1 1  In  grows  faster  than  that  of  I  In,  as  shown  in  Figure  17.  Notice  that 
family  I  In  does  not  yield  a  smoother  scaling  function  until  the  sequence  length  reaches 
20,  and  family  1 1  In  does  not  do  so  until  the  sequence  length  reaches  18.  Figure  18 
shows  the  scaling  functions  and  their  fourth  derivatives  for  representatives  of  each  of 
these  families  corresponding  to  length-30  scaling  sequences.  The  Sobolev  exponents 
of  the  three  scaling  functions  are  4.565,  5.410,  and  6.291,  respectively.  While  the 
scaling  functions  appear  quite  similar,  the  difference  shows  up  in  the  graph  of  the 
fourth  derivative. 

Similarly,  when  M  =  4  we  examined  two  PPZ  wavelet  families,  the  maximal 
vanishing  moment  family  with  filter  length  AN 

In  ;  (tt,  N;  7r/2,  N-,  37r/2,  N)  (generalized  Daubechies), 
and  a  new  family  with  filter  length  4N  —  1 

IIn  :  (tt,  iV  -  1;  ^2,  N  -  1-  37r/2,  W  -  1;  47r/5,  1;  67r/5,  1)  . 

The  new  PPZ  family  has  one  vanishing  moment  relaxed,  with  the  resulting  degrees 
of  freedom  used  to  assert  a  zero  at  the  periodic  point  47r/5.  The  scaling  functions 
from  this  family  proved  to  have  Sobolev  regularity  that  grew  significantly  faster  than 
for  the  maximal  vanishing  moment  family  In,  as  shown  in  Figure  19. 

For  M  =  3  we  define  the  families: 

In  •  (27r/3,  N]  47r/3,  N)  (generalized  Daubechies,  sequence  length  dN) 

IIn  '  (27r/3,  N  —  V,  47r/3,  iV  —  1;  tt,  1)  (sequence  length  ZN  —  1) 

and  find  that  the  asymptotic  smoothness  of  family  In  is  logarithmic  in  N  (as  men¬ 
tioned  above),  and  algebraic  in  N  (of  order  N'^^)  for  family  IIn-  Figure  20  shows  the 
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Three  PPZ  families  of  M=2  orthonormal  wavelets 

I 

0  =  Daubechies  maximal  v.m. 

+  =  Additional  zero  at  2*pi/3 

X 

X  =  Additional  zeros  at  2*pi/3  and  5*pi/6 

X  4. 

X  + 

X  + 

X  + 

X  + 

X  + 

X  +  o  ^ 

X  +  n  O 


X  +  o  ^ 

X  +  o 
6^ 


o  S 
O  ^ 


o  +  X 

o  +  V 


10 


20  30  40 

Sequence  length 


50 


60 


Figure  17:  Sobolev  exponents  for  the  PPZ  families  7/yv,  and  III^  of  M 
orthonormal  wavelets. 


=  2 


Sobolev  exponents  for  scaling  functions  from  the  new  family  and  the  minimal-length 
family.  Scaling  functions  from  the  two  M  ==  3  families  with  similar  support  lengths 

(17  and  17.5)  are  compared  in  Figure  21. 

These  examples  illustrate  the  relation  between  smoothness  and  the  specification 
of  the  Fourier  transform  of  the  filters  at  specified  preperiodic  points,  as  pointed  out 
by  previous  authors  [20],  [24].  The  patterns  described  indicate  that  a  family  of  PPZ 
wavelets  with  vanishing  specified  at  preperiodic  points  other  than  just  the  preimages 
of  27r  under  t  may,  in  general,  lead  to  greater  smoothness  than  the  maximal  vanishing 

moment  family. 
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Minimal  length  30-coeff  s.f. 


0  10  20 

Two  addl.  zeros  30-coeff  s.f. 


Fourth  derivative  of  minlength  s.f. 


Fourth  derivative  of  addl.  zero  s.f. 


Fourth  derivative  of  two  addl.  zero  s.f. 


Figure  18:  Dyadic  (M  =  2)  scaling  functions  and  their  fourth  derivatives  for  30- 
coefiicient  representatives  of  the  families  In,  Hni  S'lid  1 1  In-  From  top  to  bottom: 
Daubechies  minimum  length,  single  additional  zero  at  27r /3,  and  two  additional  zeros 
at  27r/3  and  Stt/O. 
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Figure  19:  Sobolev  exponents  for  the  families  In  and  IIn  of  M  =  4  wavelets. 


Two  PPZ  families  of  M-3  orthonormal  wavelets 


Sequence  length 


Figure  20;  Sobolev  exponents  for  the  two  families  In  and  I  In  of  M  -  3  scaling 
functions,  iV  =  1,  2,  . . . ,  15. 


M=3  N=12  minimal  length  scaling  function;  L=36,  s  =  1 .5741 


M=3  N=11  addl.  zero  at  pi  scaling  function;  L=35,  s  =  3.9383 


Figure  21:  Two  M  —  3  scaling  functions,  one  with  support  18  and  maximal  approx¬ 
imation  degree  (12),  and  the  other  with  support  17|,  approximation  degree  11,  and 
an  additional  zero  at  tt  for  the  symbol. 
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4.4.7  Quadratic-Constrained  Optimization  for  Wavelet  Filter  Design 

In  addition  to  developing  the  approaches  to  wavelet  filter  design  described  above, 
we  have  implemented  a  high-performance  design  approach  due  to  Professor  Truong 
Q.  Nguyen  of  the  University  of  Wisconsin,  a  consultant  to  Aware  on  this  contract. 
Nguyen’s  QCLS  (quadratic-constrained  least-squares)  method  [69],  [70]  is  applicable 
to  a  broad  class  of  filter  design  problems,  including  a  wide  variety  of  wavelet  fil¬ 
ter  banks.  This  approach  formulates  the  orthogonality  condition  (2)  as  a  quadratic 
constraint  on  the  wavelet  filter  coefficients,  and  optimizes  a  functional  such  as  out-of- 
band  energy  over  the  coefficient  space  (i.e.  using  the  time-domain  representation  of 
the  filter).  It  is  possible  to  include  vanishing  wavelet  moments  as  additional  quadratic 
constraints.  QCLS  has  produced  the  highest-performance  (best  stopband  attenuation 
or  subchannelization)  design  examples  known  for  cosine-modulated  filter  banks  (co¬ 
sine  packets).  We  describe  this  method  briefly  now. 

The  QCLS  algorithm  begins  with  the  vector  h  of  unknown  lowpass  filter  coef¬ 
ficients,  and  formulates  the  perfect-reconstruction  or  orthogonality  conditions  as  a 
quadratic  form  in  h.  Depending  on  the  structure  of  the  rank  M  wavelet  matrix  (arbi¬ 
trary  paraunitary  matrix,  linear-phase,  or  cosine-modulation,  as  described  m  section 
4.5  below),  the  number  Nc  of  constraints  and  their  particular  form  will  vary.  In  all 
cases,  the  Nc  independent  PR  constraints  in  (2)  can  be  expressed  as 

h‘ h  =  f  =  1,  2,  . . . ,  yVe  . 

The  objective  function  $  for  the  optimization  (taken  to  be  the  stopband  error  of 
H{z))  can  also  be  expressed  as  a  quadratic  form  in  h  via  the  eigenfilter  formulation 

[114]: 

$(h)=  r  \H{e^‘^)\^du  =  h^Ph  . 

JoJs 

Both  the  objective  function  and  constraints  admit  closed-form  expressions  for  the  gra¬ 
dient  and  Hessian,  enabling  efficient  solution  of  the  constrained  minimization  problem 

minh  $(h)  subject  to  h*  h  =  ,  1  <  ^  .  (35) 

A  nonlinear  constrained-optimization  algorithm  such  as  that  of  Schittkowski  [91] 
is  used  to  solve  (35).  The  solution  h  is  the  time-domain  impulse  response  of  a  low- 
pass  filter  with  minimal  stopband  energy,  satisfying  the  appropriate  rank  M  wavelet 

orthogonality  conditions.  j  i  i.  j 

We  have  implemented  the  QCLS  approach  for  the  design  of  cosine-modulated 
rank  M  wavelet  matrices  in  our  WaveTool  software,  as  well  as  the  earlier  lattice- 
parametrization  method  which  uses  an  unconstrained  optimization  over  a  highly  non¬ 
linear  parameter  space.  The  QCLS  method  systematically  yields  higher-performance 
filters,  as  shown  in  the  two  WaveTool  plots  below  (Figures  22  and  23). 
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Figure  22:  Frequency  response  of  an  M  =  8  length  64  lowpass  filter  for  cosine  mod 
ulation,  designed  with  QCLS  using  WaveTool  software. 


Figure  23:  Frequency  response  of  an  M  =  8  length  64  lowpass  filter  for  cosine  mod¬ 
ulation,  designed  with  lattice  parametrization  using  WaveTool  software. 
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4.5  Construction  of  Full  Rank  M  Wavelet  Matrices 

In  the  previous  section  we  saw  a  number  of  methods  for  constructing  rank  M  scaling 
sequences  of  degree  N  and  arbitrary  length.  We  now  turn  to  the  design  of  the 
corresponding  wavelet  sequences,  i.e.  to  the  problem  of  constructing  a  complete 
rank  M  wavelet  matrix  given  its  first  row.  In  the  rank  M  case,  there  is  considerable 
freedom  in  such  a  construction.  We  describe  here  three  distinct  methods.  The  first 
is  completely  general,  and  is  based  on  the  prior  choice  of  a  wavelet  lowpass  filter 
(scaling  sequence)  and  a  Haar  wavelet  matrix  [46].  The  second  and  third  methods 
involve  more  specialized  structures  -  linear  phase  rank  M  wavelet  matrices  and  rank 
M  cosine-modulated  wavelets,  respectively.  All  three  approaches  have  associated  fast 
algorithms  for  computation. 

4.5.1  General  Haar-based  constructions 

First,  we  present  a  method  [39]  for  constructing  a  full  wavelet  matrix  given  its  first 
row  (the  scaling  sequence)  and  an  Af  x  M  matrix  which  we  call  the  characteristic 
Haar  matrix  [46]. 

A  Haar  wavelet  matrix  is  an  orthogonal  matrix  (up  to  scalar  multiplication)  whose 
first  row  is  all  ones;  that  is,  its  entries  hs^k  satisfy 

k 

and  ho,k  =  1  V/: . 

Observe  that  every  such  Haar  wavelet  matrix  is  a  rank  M  wavelet  matrix  (with 
K  =  M)  and  that  every  M  x  M  wavelet  matrix  is  a  Haar  matrix.  Useful  examples  of 
Haar  wavelet  matrices  include  the  M-point  FFT,  type  H  Discrete  Cosine  Transform, 
and  Hadamard  matrix.  The  collection  of  rank  M  Haar  matrices  is  isomorphic  to  the 
group  of  orthogonal  matrices  of  rank  M  —  1. 

It  will  serve  us  to  think  of  our  wavelet  matrices  as  being  M  x  Mg  for  some  integer 
g',  we  can  always  pad  each  row  with  zeros  to  bring  the  wavelet  matrix  into  this  form. 
g  is  called  the  genus  of  the  wavelet  matrix.  If  we  break  up  the  M  x  Mg  wavelet 
matrix  A  into  its  constituent  M  x  M  blocks 

A  =  (Ao  Ai  ...  A,^i)  (36) 

then  the  characteristic  Haar  matrix  associated  with  A  is  given  by 

Ho  =  Ao  +  Ai  +  A^_i  . 

It  can  be  checked  that  Hq  is  in  fact  a  Haar  wavelet  matrix. 

In  [46]  we  solved  the  following  problem:  given  a  Haar  wavelet  matrix  Ho  and 
a  scaling  sequence  ao,  construct  a  full  wavelet  matrix  A  whose  first  row  is  ao  and 
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whose  characteristic  Haar  matrix  is  Hq.  That  explicit  construction  can  be  clarified 
and  refined  using  Vaidyanathan’s  paraunitary  factorization  technique.  In  particular, 
a  parametrization  of  the  choice  of  the  M  —  1  wavelets  is  given  by  the  choice  of  the 
characteristic  Haar  matrix. 

Working  in  the  ^-transform  domain,  Vaidyanathan  [116]  has  proven  that  every 
paraunitary  polyphase  matrix  H(0)  of  McMillan  degree^  K  can  be  factored  into  the 
form 

=  f  n(I  -  Ho  ,  (37) 

\Ai=0  / 

where  each  is  a  unit  M-vector.  H(z)  will  be  the  polyphase  matrix  of  a  wavelet 
matrix  if  and  only  if  Ho  is  a  Haar  wavelet  matrix,  so  (37)  provides  a  factorization  of 
all  wavelet  matrices  with  polyphase  matrix  of  McMillan  degree  K.  We  refer  to  the 
term 

I  -  VfcV^  +  ZWkvl 

as  a  prime  factor  of  the  polyphase  matrix. 


Theorem  4.6  Given  a  scaling  sequence  ao  of  overlap  g  and  a  characteristic  Haar 
matrixUo,  there  exists  a  unique  wavelet  matrix  of  McMillan  degree  g  —  1  {i.e.  with  g— 
1  prime  factors)  whose  first  row  is  ao  and  with  characteristic  Haar  Ho-  Furthermore, 
we  can  explicitly  construct  the  prime  factors  I  —  va;V](  +  zvkvl  {and  thus  the  wavelet 
matrix)  from  ao  and  Hq. 


The  proof  of  this  theorem  (including  explicit  construction  of  the  prime  factors) 
appears  in  [39]. 

As  an  example,  we  use  this  method  to  construct  a  wavelet  matrix  with  M  =  4 
and  g  =  2,  with  approximation  degree  N  =  g  —  2.  The  minimal  length  lowpass  filter 
(scaling  sequence)  for  this  case  was  given  previously  (equation  (20)).  The  full  wavelet 
matrix  constructed  using  Theorem  4.6,  with  the  sequence  (20)for  its  first  row  and  the 
rank-4  DCT-II  for  its  characteristic  Haar  matrix  Hq  is 
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polyphase  matrix  of  McMillan  degree  K  will  correspond  to  a  wavelet  matrix  of  overlap 
while  a  wavelet  matrix  of  overlap  K  1  has  a  polyphase  matrix  with  McMillan  degree  at  least  K. 
However,  there  exist  wavelet  matrices  of  overlap  K  +l  and  McMillan  degree  strictly  greater  than  K; 
for  examples  see  [50] .  The  construction  summarized  here  and  presented  in  detail  in  [39]  describes  a 
unique  wavelet  matrix  with  first  row  uq  and  characteristic  Haar  Ho  and  having  a  polyphase  matrix 
of  McMillan  degree  K. 
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The  frequency  responses  of  the  lowpass  filter  (scaling  sequence)  and  the  three  wavelet 
filters  (i.e.,  the  four  rows  of  the  naatrix)  appear  in  Figure  24. 


Figure  24:  Magnitude  frequency  responses  of  filters  in  a  minimal  length  M  4  ,  N 
2  wavelet  matrix  based  on  a  DCT-II  characteristic  Haar  matrix. 

When  the  rank  M  >  2,  the  construction  of  Theorem  4.6  gives  considerable  free¬ 
dom  in  the  choice  of  the  M  -  1  wavelet  sequences.  We  have  exploited  the  range  of 
possibilities  in  this  parametrization  for  image  compression  applications,  as  reported 
in  [45]. 

4.5.2  Linear-Phase  Rank  M  Wavelets,  or  GenLOT  constructions 

A  more  restrictive  subclass  of  the  rank  M  wavelets  are  those  with  linear-phase,  i.e. 
such  that  each  of  the  M  filters  have  linear-phase  symmetry.  Such  symmetry  is  of  great 
utility  in  applications  to  image  compression,  where  symmetric  extension  of  data  is 
the  preferred  method  for  handling  image  boundaries  [95].  Linear-phase  filters  are  also 
of  value  in  a  wavelet  transform  because  they  preserve  centers  of  mass  m  an  iterated 
signal  decomposition.  In  the  rank  2  case,  wavelet  matrices  cannot  simultaneously 
satisfy  orthogonality  (2)  and  symmetry,  other  than  the  nearly  trivial  Haar  exaniple. 
This  has  been  overcome  by  designing  biorthogonal  rank  2  wavelet  matrices  [73],  [11], 
[117].  However,  when  M  >  2,  linear-phase  and  orthogonality  can  be  simultaneously 
satisfied.  Such  orthogonal  linear-phase  filter  banks  have  recently  been  parametrized 
[97],  [79],  at  least  when  M  is  even.  We  have  used  these  parametrizations  to  deter¬ 
mine  the  M-band  orthogonal  linear-phase  wavelets,  and  successfully  applied  the  new 
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wavelet  matrices  to  problems  in  image  compression. 

Let  us  briefly  summarize  the  parametrizations  of  Soman  and  deQueiroz.  One  of 
these  works  [79]  identifies  M-band  orthonormal  filter  banks  in  which  each  filter  has 
linear-phase  symmetry  as  generalizations  of  the  lapped  orthogonal  transform  (LOT) 
[66],  and  confers  the  name  GenLOT.  As  described  previously,  a  rank  M  wavelet 
matrix  can  be  described  in  the  z-transform  domain  by  its  polyphase  matrix  E(.2); 
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When  the  rank  M  is  even,  and  the  filters  hk[n]  (with  ^-transforms  Hfc(z))  have  length 
MN^  the  polyphase  matrix  has  the  form 


E{z)  =  KN-i{z)KN-2{^) . . .  Ki(z)Eo 


(38) 


where  each 


Ki(2) 


■  Ui 

0 

■  I  I 

■  I  o' 

I  I 

0 

I  -I 

0  z-n 

I  -I 

I  is  the  rank  M/2  identity  matrix,  while  U,-  and  V,-  are  arbitrary  rank  M/2  orthogonal 
matrices.  Eq  is  a  rank  M  unitary  matrix  with  M/2  symmetricand  M/2  antisymmetric 
rows  (such  as  the  DOT- IV);  it  can  be  factored  as 
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J  is  the  familiar  reverse  identity  matrix  (of  rank  M/2),  while  Do  and  Di  are  arbitrary 
orthogonal  matrices  of  rank  M/2. 

This  GenLOT  formulation  can  be  used  to  construct  of  M-band  linear-phase  or¬ 
thonormal  wavelets.  First,  all  possible  GenLOTs  with  one  vanishing  moment  may 
be  described  in  terms  of  certain  rotation  matrices.  As  described  in  Theorem  4.1,  an 
orthogonal  filter  bank  has  one  vanishing  moment  if  it  satisfies 
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for  z  —  1.  This  is  equivalent  to  the  linear  equation  (3)  in  the  definition  of  a  wavelet 
matrix.  We  find  that  a  GenLOT  with  N  —  1  factors  Eii{z)  will  have  one  vanishing 
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moment  if  and  only  if  the  rotation  matrices  U,-  satisfy 


Uyv-lUiV-2 


..UiDo 


■  1  ■ 

1 - 

CM 

1 

0 

.  1 . 

0 

Clearly,  the  number  of  free  parameters  available  for  filter  design  (to  optimize  such 
properties  as  stopband  attenuation  or  coding  gain)  increases  with  both  the  number 
of  channels  M  and  the  filter  length,  which  is  determined  by  the  number  of  factors 
N.  With  even  the  simplest  iV  =  1  LOTs,  it  is  possible  to  create  a  filter  bank  based 
on  the  DCT-IV  having  one  vanishing  moment  for  use  in  a  wavelet  decomposition.  A 
4-band  GenLOT  with  N  =  2  (filter  length  12)  and  one  vanishing  moment  is  shown 
in  Figure  25. 


Figure  25:  Magnitude  responses  (in  dB)  of  a  rank  4  genus  3  linear-phase  wavelet 
matrix  (GenLOT)  with  one  vanishing  moment. 

It  is  of  interest  to  create  wavelet  filters  with  more  than  one  vanishing  moment; 
indeed,  the  interpolation/approximation  properties  of  such  filters  lead  to  their  supe¬ 
riority  for  wavelet-based  image  coding  systems  [63].  An  rank  M  wavelet  matrix  will 
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have  two  vanishing  moments  if  it  satisfies  (39)  as  well  as  the  second-order  condition 
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at  z  =  1.  In  the  GenLOT  case,  this  leads  to  a  set  of  linear  constraint  equations  on 
the  rotation  matrices  Ui,  and  that  parametrize  the  wavelet  matrix.  These 
constraints  can  be  solved  to  reduce  the  parameter  space,  enabling  one  to  optimize  over 
the  remaining  parameters  to  design  a  rank  M  linear-phase  wavelet  matrix  with  two 
vanishing  moments  and  other  desirable  properties  such  as  high  stopband  attenuation 
or  maximal  smoothness.  The  frequency  responses  of  such  a  wavelet  matrix  with 
M  —  A  are  shown  in  Figure  26. 


Figure  26:  Magnitude  responses  (in  dB)  of  a  rank  4  genus  4  linear-phase  wavelet 
matrix  (GenLOT)  with  two  vanishing  moments  and  maximal  Sobolev  smoothness. 

We  applied  these  new  higher  rank,  linear-phase  wavelets  to  image  compression  [44], 
with  results  superior  to  those  obtained  using  rank  2  wavelets.  Particular  advantages 
were  found  in  the  compression  of  fingerprints  and  seismic  data.  Rank  4  wavelet 
matrices  appear  naturally  in  the  special  wavelet-packet  tree  (Figure  27)  specified  by 
the  FBI  for  their  Wavelet  Scalar  Quantization  fingerprint  compression  algorithm  [7]; 
this  transform  tree  can  be  obtained  as  the  cascade  of  two  rank  4  wavelet  matrices, 
followed  by  a  rank  2  wavelet  in  the  lowest-frequency  subband.  When  we  substituted 
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the  GenLOT  examples  given  above  into  the  transform  structure,  lower  maximum 
(Chebyshev)  errors  resulted  across  all  compression  ratios;  the  effect  was  particularly 
pronounced  at  high  compression  ratios  (low  bitrates). 


Figure  27:  Wavelet  transform  tree  specified  by  FBI  for  fingerprint  compression. 


4.5.3  Cosine-modulation 


One  of  the  most  effective  methods  for  the  design  and  implementation  of  rank  M  wave¬ 
let  filter  bank  transforms  is  via  cosine  modulation.  In  this  approach,  a  single  lowpass 
“prototype  filter”  is  modulated  to  create  a  complete  rank  M  wavelet  matrix,  using 
a  DCT-based  modulation  matrix.  This  method  yields  both  the  highest-performance 
wavelet  filters,  and  the  fastest  known  algorithms  for  rank  M  wavelet  transform  compu¬ 
tation.  We  summarize  the  cosine  modulation  construction  here;  such  wavelet  matrices 
also  go  by  the  name  of  cosine  packets  and  discrete  local  cosine  transforms. 

In  the  most  general  biorthogonal  case  with  DOT- IV  modulation,  the  impulse 
responses  of  the  analysis  filters  hk[n]  are  cosine-modulated  versions  of  a  prototype 
filter  h[n]  of  length  Nk,  and  the  synthesis  filters  can  also  be  obtained  via  cosine- 
modulation  of  a  length  Nf  prototype  filter  /[n].  The  overall  reconstruction  delay  D 
of  the  filter  bank  can  be  fixed  arbitrarily  in  the  range  D  G  [0,  N/  +  Nh  -  1]-  For  a 
given  delay  D  =  2sM  +  d  (where  0  <  d  <  2M),  the  relation  between  the  analysis 
filters,  the  synthesis  filters,  and  their  prototypes  can  be  stated  as  follows: 


hk[n]  = 
fk[n]  = 


2h[n\  cos 
2/[n]  cos 


(2fc  + 1)^(1 -|)  +  «i 
(2k  +  l)^(n  -  |)  -  h 


(41) 

(42) 


with  Ok  =  (-l)^f-  Note  that  the  modulation  does  not  depend  on  the  filter  length 
but  only  on  the  delay  of  the  system.  However,  if  we  restrict  ourselves  to  the  case 
where  h[n]  =  /[n]  and  both  are  length  N  linear-phase  prototype  filters,  the  delay  is 
constrained  to  be  D  =  N  —  1  and  this  is  the  same  modulation  as  in  [57]. 
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As  an  example,  the  lowpass  prototype  filter  shown  in  Figure  22  can  be  modulated 
to  yield  the  rank  8  wavelet  matrix  (filter  bank)  shown  in  Figure  28  below.  Notice  that 
the  superior  stopband  attenuation  of  the  prototype  filter  is  preserved  in  the  modulated 
wavelet  matrix.  The  superior  sub  channelization  properties  of  cosine- modulation  is 
due  to  a  combination  of  the  relatively  clean  parametrization  of  the  class  of  possible 
prototype  filters  and  the  high-performance  QCLS  algorithm. 


Figure  28:  Frequency  responses  of  an  M  =  8  length  64  cosine-modulated  wavelet 
matrix,  designed  using  WaveTool  software. 

Further  details  on  cosine-modulated  filter  banks  and  wavelets  appear  in  a  number 
of  references,  including  [66],  [57],  [81],  [102],  [36]  [42]  among  others. 
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5  WaveTool  Software 


Foremost  among  our  accomplishments  under  this  contract  was  the  development  and 
release  of  the  WaveTool  software.  This  is  a  UNIX-based  software  tool  for  rapidly 
prototyping  and  evaluating  general  higher-rank  wavelet  decompositions,  with  a  so¬ 
phisticated  graphical  user  interface.  The  WaveTool  grew  out  of  our  early  experience 
of  applying  wavelet  methods  here  at  Aware.  We  used  to  laboriously  hand-code  rou¬ 
tines  to  perform  wavelet/multirate  decompositions  for  each  new  application  as  it 
arose.  We  had  one  set  of  software  which  did  rank  2  wavelet  decompositions  for  image 
compression,  another  which  performed  a  rank  32  modulated  lapped  transform  for 
CD-quality  audio  compression,  and  a  third  for  transient  detection.  While  each  new 
application  might  require  a  new  basis  (new  filters  or  new  tree),  the  underlying  me¬ 
chanics  of  the  decomposition  was  always  the  same.  Furthermore,  a  significant  step  in 
the  development  of  wavelet  algorithms  proved  to  be  the  hunt  for  the  right  basis.  Out 
of  this  experience  grew  our  specification  for  the  WaveTool  -  a  piece  of  software  which 
would  enable  the  user  to  rapidly  prototype  wavelet  algorithms  with  a  wide  range  of 
filter  and  tree  choices,  and  underlying  fast  computer  programs  for  implementing  the 
wavelet  decompositions,  once  chosen. 

The  WaveTool  software  provides  the  capability  to  design  a  variety  of  wavelet  and 
multirate  filter  banks,  assemble  them  into  arbitrary  tree  structures,  and  process  data 
through  the  system  in  either  subband  (one-signal-into-many-subbands)  or  transmul¬ 
tiplexor  (many-signals-into-one)  modes.  The  WaveTool  also  offers  extensive  graphical 
display  capabilities  for  both  filter  bank  components  and  input  and  output  data.  Per¬ 
haps  the  easiest  way  to  describe  the  software  is  to  take  the  reader  on  a  tour  of  its 
functionality. 

The  heart  of  the  WaveTool  is  its  Tree  Window,  in  which  the  user  interactively 
defines  a  tree  structure  of  paraunitary  filter  banks.  He  or  she  does  so  by  picking 
off  a  menu;  the  user  may  pick  from  a  library  of  predesigned  filter  banks,  may  read 
in  a  filter  bank  of  his  own  design,  or  may  interactively  design  a  filter  bank  with  the 
software.  There  are  choices  for  either  wavelet  or  Modulated  Lapped  Transform  (MLT) 
designs.  In  the  case  of  a  wavelet  filter  bank,  the  user  specifies  two  parameters:  the 
rank  or  number  of  channels  M  and  the  filter  length  L.  The  resulting  orthonormal 
(paraunitary)  wavelet  filter  has  the  maximal  number  of  vanishing  wavelet  moments 
for  the  given  filter  length;  these  filter  banks  generalize  Daubechies  construction  to 
the  M-channel  case,  as  described  in  section  4.4.1.  For  M  >  2,  the  wavelet  filters  are 
determined  by  the  characteristic  Haar  matrix  (polyphase  matrix  evaluated  at  z  =  1); 
this  is  set  to  the  discrete  cosine  transform  (DCT),  because  this  choice  has  proven 
effective  in  applications.  The  other  possibility  for  designing  a  paraunitary  filter  bank 
is  the  cosine-modulated  or  MLT  case,  in  which  one  designs  a  prototype  lowpass  filter 
which  is  modulated  via  a  DCT  to  uniformly  partition  the  full  frequency  spectrum 
(section  4.5.3).  Again,  the  user  prescribes  the  rank  M  and  the  overall  filter  length 
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Figure  29:  Two  MLT  prototype  filter  designs  in  WaveTool,  with  M  =  8,  L  =  32. 


L.  When  the  filter  length  is  greater  than  2M,  there  are  free  parameters  which  may 
be  used  to  optimize  the  filter  for  superior  stopband  attenuation,  narrow  transition 
bandwidth,  etc.  Figure  29  shows  two  MLT  prototype  filter  designs  from  the  WaveTool 
for  rank  M  =  8  and  filter  length  L  =  82  which  trade  off  between  height  of  the  first 
sidelobe  and  transition  bandwidth. 

Once  the  user  has  designed  a  filter,  he  may  then  “accept”  it  into  the  nascent  tree 
structure  in  the  Tree  Window.  The  filter  bank  coefficients  (time  domain  impulse 
responses)  are  also  saved  to  an  ASCII  file.  The  impulse  and  frequency  responses  of 
any  filter  bank  in  the  tree  may  be  displayed  with  a  mouse  click. 

A  salient  feature  of  the  WaveTool  is  the  ability  to  create  arbitrary  tree  structures 
of  paraunitary  filter  banks.  The  classical  wavelet /multiresolution  decomposition  [23], 
[64]  is  obtained  by  cascading  rank  2  filter  banks  on  their  lowpass  outputs  only.  Given 
any  rank  M  filter  bank  in  the  Tree  Window,  the  user  may  create  such  a  Mallat  tree 
with  a  single  menu  pick.  However,  non-Mallat  trees  have  also  proved  quite  useful  in 
signal  analysis  (e.g.  wavelet  packet  decompositions  [22],  the  FBI’s  WSQ  fingerprint 
compression  standard  [7],  and  audio  compression  using  a  psychoacoustic  model  [98]). 
Any  such  tree  may  be  created  with  the  WaveTool. 

Once  the  user  has  constructed  a  tree-structured  filter  bank,  he  or  she  can  use  it 
to  operate  on  data.  One  loads  in  an  input  data  file  (in  ASCII,  binary,  or  MATLAB 
format);  the  input  signal  may  be  displayed  in  a  graphics  window.  For  example,  we 
have  loaded  in  the  piecewise  quadratic  signal  “nuquad”  and  displayed  it,  as  shown  in 
Figure  30.  Plots  such  as  this  may  be  zoomed  and  relabeled  at  will,  and  dumped  to 
PostScript  output  for  inclusion  in  reports. 
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Figure  30:  Piecewise  quadratic  input  signal. 

We  then  pick  an  “M2L6”  wavelet  filter  from  the  Pick  Filter  menu  (i.e.  Daube-chies’ 
rank  2,  6-tap  orthonormal  wavelet  filter  bank),  and  cascade  it  into  a  three-level  Mallat 
tree,  as  shown  in  Figure  31.  Notice  that  the  Tree  Window  uses  the  convention  that 
lowpass  outputs  are  toward  the  top  of  the  window.  When  we  hit  the  Analyze  button, 
we  perform  a  decomposition  of  the  input  into  the  octave-based  wavelet/subband 
decomposition.  We  can  display  the  subband  outputs  as  streaming  from  the  output 
leaves  of  the  tree-structured  filter  bank  by  hitting  the  Show  Subbands  button;  this  is 
also  shown  in  Figure  31.  The  polynomial  interpolation  properties  of  the  Daubechies 
6-tap  filter  are  shown  clearly  here  -  the  only  nonzero  bandpass/highpass  subband 
outputs  occur  at  the  knots  of  the  piecewise  quadratic  input. 
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nuquad 


Figure  31:  “D6”  wavelet  decomposition  of  the  piecewise  quadratic  signal. 
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Figure  32:  Piecewise  quadratic  input  signal,  and  its  coarse  reconstruction. 

Having  applied  the  wavelet  transform,  one  can  then  output  the  transform  coeffi¬ 
cients  to  ASCII  or  MATLAB  data  files.  The  MATLAB  output  option  is  particularly 
useful,  for  one  can  load  the  transform  coefficients  in  as  matrix  variables  and  perform 
quantization,  peak-picking,  and  other  operations  on  the  coefficients  before  saving  back 
to  file,  and  reloading  into  the  WaveTool  for  the  synthesis  (reconstruction)  operation. 
For  example,  let  us  take  the  wavelet  subbands  from  our  D6-Mallat  tree  decomposi¬ 
tion  of  the  piecewise  quadratic  signal,  load  them  into  MATLAB,  and  zero  out  the 
bandpass  and  highpass  subbands.  We  then  load  the  remaining  coarse-scale  transform 
coefficients  back  into  the  WaveTool  and  reconstruct  into  a  signal.  Figure  32  shows  the 
original  signal,  followed  by  the  reconstruction  from  the  coarse-scale  approximation  at 
a  delay  of  28  samples.  The  coarse  reconstruction  hardly  differs  from  the  original  piece- 
wise  quadratic  input,  again  because  of  the  interpolation  properties  of  the  Daubechies 
filters.  Furthermore,  one  may  measure  signal  reconstruction  error  in  the  WaveTool, 
using  a  variety  of  norms  For  example,  the  coarse  reconstruction  of  Figure 

32  has  an  £°°  error  of  0.1264093.  These  data  input/output  features  are  not  the  only 
means  by  which  the  WaveTool  interacts  with  other  software;  once  a  wavelet/subband 
decomposition  has  been  designed  in  the  WaveTool,  the  algorithm  may  be  incorpo¬ 
rated  into  a  larger  MATLAB  simulation  (of,  say,  a  compression  or  communications 
system)  via  MEX-file  implementations  of  the  filtering  algorithms.  Further  details  on 
the  operation  of  the  WaveTool  may  be  found  in  the  software’s  User’s  Guide  [4]  and 
Reference  Manual  [3]. 

In  addition  to  its  use  in-house  at  Aware  for  tasks  such  as  waveform  design  for 
Discrete  Wavelet  MultiTone  modulation  (section  6.1),  the  WaveTool  software  has 
been  installed  at  a  number  of  academic,  government,  and  industrial  beta  sites,  and 
ultimately  released  as  a  commercial  software  product.  This  is  discussed  in  section 
10.1  below. 
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6  Applications 

We  now  turn  our  discussion  to  the  third  leg  of  this  project  -  application  of  higher-rank 
wavelet  techniques  to  telecommunications  and  signal  processing.  A  surprise  winner 
has  been  the  use  of  rank  M  wavelet  filter  banks  to  multicarrier  modulation,  a  leading 
technique  for  high-bitrate  last-mile  telecommunications,  which  we  discuss  first.  We 
then  examine  the  application  of  the  wavelet  transforms  and  design  techniques  to  the 
compression  of  image  data,  including  new  data  regimes  such  as  seismic,  multispec- 
tral,  and  sonar  data.  Finally,  we  review  our  application  of  emerging  multiwavelet 
techniques  to  signal  and  image  processing. 

6.1  DWMT  -  Multicarrier  Modulation  via  Rank  M  Wavelets 

Multicarrier  modulation  has  recently  emerged  as  a  superior  method  for  achieving  high- 
bitrate  transmission  over  the  “last  mile”  of  the  installed  wire  plant.  This  application 
is  of  particular  importance,  given  the  dramatic  increase  in  use  of  the  Internet  and  the 
attendant  need  for  high-bandwidth  connections  to  the  home.  The  last  mile  usually 
takes  the  form  of  either  twisted-pair  copper  wire  or  a  hybrid  fiber-coax  cable.  ANSI’s 
T1E1.4  committee  has  chosen  a  multicarrier  scheme  as  the  standard  for  Asymmetric 
Digital  Subscriber  Line  (ADSL)  transmission.  ADSL  is  an  emerging  technology  that 
promises  upwards  of  6  Mbits/sec  into  the  home  or  business,  using  the  existing  copper 
plant,  while  simultaneously  allowing  conventional  telephone  service  over  the  same 
line.  Multicarrier  can  also  be  used  for  hybrid  fiber-coax  networks,  to  provide  two-way 
communications  services  (data  delivery  or  telephony)  over  the  existing  cable  television 
(CATV)  network. 

Bitrates  in  the  Megabit/sec  range  employ  frequency  bandwidths  of  1  MHz  or  more; 
both  twisted-pair  and  HFC  have  particular  impairments  across  these  wide  frequency 
bands.  In  the  case  of  twisted-pair  copper,  the  available  signal-to-noise  ratio  (SNR) 
for  information  throughput  varies  strongly  as  a  function  of  frequency,  particularly 
over  long  (10,000  feet  or  more)  loops  and  those  with  “bridge  taps”  (unterminated, 
unused  stubs  of  wire).  Both  twisted  pair  and  HFC,  but  particularly  the  HFC  up¬ 
stream  channel  (5-40  MHz)  are  also  susceptible  to  powerful  narrowband  interference, 
such  as  that  caused  by  AM  radio  transmissions.  As  we  discuss  below,  each  of  these 
impairments  can  be  bypassed  or  overcome  by  a  multicarrier  modulation  scheme. 

In  multicarrier  modulation,  the  transmission  channel  is  partitioned  into  a  number 
M  of  subchannels  (usually  128  <  M  <  512),  each  with  its  own  associated  carrier 
[78].  This  is  usually  accomplished  digitally,  using  an  orthogonal  transformation.  At 
the  receiver,  the  inverse  transform  is  performed  to  demodulate  the  data.  This  can 
be  interpreted  as  transmultiplexer  taking  time-division-multiplexed  (TDM)  data  into 
frequency-division-multiplexed  (FDM)  data.  Thus  multicarrier  modulation  provides 
an  efficient  means  to  access,  transmit,  and  distribute  multiple  data  streams.  Each 
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subchannel  has  its  own  characteristic  SNR  that  can  be  measured  and  is  then  used 
to  determine  the  constellation  size  (number  of  bits  supported)  for  that  subchannel. 
In  contrast  to  a  single-carrier  scheme,  multicarrier  gives  a  fine-grained  (M  grams,  to 
be  exact)  decomposition  of  the  transmission  channel  for  the  purpose  of  assigning  bit 
levels.  This  allows  for  more  nearly  optimal  use  of  the  available  bandwidth  in  the  case 
of  SNR  that  varies  strongly  with  frequency,  as  in  twisted-pair  transmission.  Multi¬ 
carrier  also  provides  superior  immunity  to  impulse  noise  (compared  to  single-carrier 
systems),  and  is  also  particularly  elTective  at  combatting  narrowband  interference.  A 
subchannel  that  is  affected  by  a  narrowband  interferer  can  be  easily  identified,  and 
the  data  rate  in  that  subchannel  reduced. 

Early  multicarrier  schemes  [124],  [54]  used  a  Fourier  transform  (DFT)  as  the 
orthogonal  modulating  transform.  This  has  the  advantage  of  fast  computational 
algorithms,  but  also  the  weakness  of  large  spectral  overlap  among  the  frequency 
responses  Ak{u)  of  the  distinct  signalling  waveforms.  This  can  lead  to  substantial 
intersymbol  interference,  and  a  lack  of  robustness  to  narrowband  interference. 

Aware,  Inc.  has  pioneered  the  use  of  wavelets  and  multirate  filter  banks  for  multi 
carrier  modulation  [111],  [90]  [88].  This  approach,  called  Discrete  Wavelet  MultiTone 
(DWMT)  or  overlapped  multitone  modulation,  yields  superior  bandwidth  utiliza¬ 
tion  and  robustness  in  the  face  of  impulsive  and  narrowband  noise  [127].  The  block 
Fourier  transform  is  replaced  by  an  orthogonal  overlapped  transform  (a  rank  M  wave¬ 
let  matrix  with  genus  g  >  1).  The  wavelet  orthogonality  conditions  (2)  prove  to  be 
equivalent  to  zero  intersymbol  and  inter-channel  interference  (ISI  and  ICI)  among 
the  signalling  waveforms.  The  use  of  an  overlapped  transform  enables  a  tradeoff  be¬ 
tween  time  duration  of  the  transform  and  spectral  isolation  of  the  subchannels;  even  a 
low  degree  of  overlap  yields  significant  improvements  over  the  block  (nonoverlapping) 
Fourier  case.  The  superior  subchannelization  offered  by  the  filter  bank  leads  to  dra¬ 
matically  superior  performance  in  the  presence  of  narrowband  (e.g.  radio-frequency) 
interference. 

6.1.1  Optimization  of  wavelet  algorithms  for  DWMT 

The  various  techniques  discussed  in  section  4  for  wavelet  filter  design  and  computation 
have  a  direct  application  to  DWMT  modulation.  In  particular,  we  have  found  the 
cosine-modulated  filter  banks  (Section  4.5.3)  most  useful  as  an  overlapped  orthogonal 
transform  for  multitone  modulation.  Cosine-modulated  filter  banks  are  determined 
by  a  single  lowpass  “prototype”  filter.  This  simplifies  the  filter  bank  design  problem, 
enabling  rapid  optimization  of  the  wavelet  filter  bank  for  properties  such  as  sub¬ 
channel  isolation.  The  notion  of  genus  or  overlap  provides  a  handy  parameter  for  the 
tradeoff  between  filter  bank  latency  (and  computational  complexity)  and  narrowband 
interference  rejection,  for  a  given  number  of  subchannels  M .  This  is  illustrated  in 
Figure  33,  which  compares  two  different  wavelet  filter  bank  responses  (sets  of  tones) 


56 


CQ 

"O 


Figure  33:  Subchannel  frequency  responses  for  a  multitone  system  with  M  =  8.  Top: 
nonoverlapped  system  based  on  the  DFT.  Middle:  overlapped  multitone  system  with 
genus  g  =  4.  Bottom:  overlapped  multitone  system  with  genus  g  =  8. 

with  those  of  a  Fourier  transform. 

The  top  plot  in  Figure  33  shows  the  frequency  responses  of  the  DFT  (used  in 
conventional  discrete  multitone  modulation)  with  M  —  8.  The  peak  sidelobe  is  at 
-13  dB,  with  1//  decay.  The  second  plot  shows  the  subchannel  frequency  responses 
of  an  M  =  8  genus  g  =  4  cosine-modulated  wavelet  filter  bank,  as  used  in  DWMT. 
This  transform  has  sidelobes  below  -35  dB.  The  third  plot,  of  an  M  =  8  genus 
g  =  8  transform  shows  sidelobes  below  -50  dB,  clearly  demonstrating  the  superior 
subchannelization  offered  by  the  overlapping  wavelet  transform.  The  wavelet  filter 
banks  displayed  in  this  example  were  designed  using  the  QCLS  algorithm  (section 
4.4.7)  as  implemented  in  the  WaveTool  software. 

The  structures  of  section  4  also  lead  directly  to  fast  algorithms  for  wavelet-based 
modulation.  In  particular,  the  cosine-modulated  filter  banks  can  be  realized  as  a 
combination  of  a  Discrete  Cosine  Transform  (DCT)  and  a  “polyphase  windowing” 
step.  Making  use  of  previous  work  on  fast  algorithms  for  the  DCT,  the  entire  cosine- 
modulated  filter  bank  computation  may  be  done  at  a  cost  not  much  greater  than  that 
of  the  DFT  with  the  same  number  of  subchannels.  Details  appear  in  [110],  [65],  [66]. 
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6.1.2  Hardware  demonstration  of  a  DWMT  modem 

We  have  developed  a  prototype  DWMT  modem  and  made  laboratory  measurements 
that  demonstrate  the  superior  performance  of  the  system  in  the  presence  of  narrow- 
band  interference.  Figure  34  shows  a  block  diagram  of  the  prototype  system.  The 
input  data  consists  of  a  serial  TDM  data  stream  that  is  divided  into  frames  of  length 
B  bits.  Given  an  input  data  rate  of  Rb  and  a  multicarrier  frame  duration  T,  the 
number  of  bits  per  frame  is  .B  =  RbT.  For  the  prototype  system,  the  parameters  are 
Rb  =  2.048  Mbits/sec,  T  =  125 /^sec,  and  B  =  256  bits/frame.  Nominal  data  rates 
on  the  order  of  5  Mbits/sec  are  achievable  and  a  maximum  rate  of  approximately  8 
Mbits/sec  can  be  obtained  in  a  2  MHz  band. 

The  data  bits  are  encoded  into  multilevel  PAM  symbols  and  the  PAM  symbols 
are  then  orthogonally  mapped  to  individual  frequency  subchannels  via  the  inverse 
wavelet  transform  modulator.  Note  that  the  subchannels  are  grouped  into  pairs  and 
modulated  with  the  equivalent  of  a  QAM  constellation.  The  wavelet  transform  is 
based  on  a  single  rank  M  wavelet  matrix  (corresponding  to  M  subchannels)  of  genus 
g.  The  properties  of  the  transform  are  such  that  the  subchannels  overlap  spectrally 
and  the  pulses  transmitted  in  a  subchannel  overlap  in  time,  while  orthogonality  among 
all  symbols  is  maintained.  The  frequency  overlap  of  adjacent  subchannels  results  in 
spectrally  efficient  transmission  and  the  time  overlap  of  the  pulses  provides  spectral 
shaping  of  the  individual  subchannel  filters.  The  transform  modulator  produces  a 
time  domain  sequence  that  is  output  to  the  D/A  converter.  In  the  downstream 
direction,  the  signal  is  broadcast  to  all  modems  and  in  the  upstream  direction,  the 
analog  signals  from  the  subscriber  modems  are  power  combined  and  then  transmitted 
to  the  head  end.  The  experimental  setup  consisted  of  an  head  end  transceiver  (as 
would  be  found  at  an  Optical  Network  Unit)  and  3  remote  transceivers  (Figure  35). 

At  the  receiver,  the  baseband  analog  signal  is  filtered  and  digitized.  The  digital 
time  domain  signal  is  transformed  back  into  the  PAM  symbols  via  the  wavelet  trans¬ 
form.  The  demodulating  transform  is  the  “analysis”  filter  bank  based  on  the  same 
rank  M  genus  g  wavelet  matrix  as  the  modulating  transform.  At  the  receiver,  this 
provides  a  set  of  matched  filters  with  respect  to  the  modulating  transform.  An  equal¬ 
izer  is  used  to  correct  for  channel  distortion  and  timing  errors.  After  equalization, 
the  PAM  symbols  are  decoded  into  bits  and  the  original  data  stream  is  reproduced 
by  the  parallel-to-serial  multiplexer.  The  receiver  can  be  constructed  to  demodulate 
all  of  the  subchannels  to  recover  the  complete  data  stream  or  a  specific  subset  of 
the  subchannels  for  information  addressed  to  individual  users.  Dynamic  allocation  of 
subchannels  and  transmission  capacity  can  be  implemented  based  on  user  demand. 

The  measured  results  are  presented  in  Figures  36  through  41.  In  Figure  36,  the 
spectrum  of  the  transmitted  signal  with  64  QAM  modulation  on  the  payload  tones 
is  given.  In  addition  to  the  group  of  tones  used  for  payload  (P),  a  set  of  tones  at  the 
upper  and  lower  band  edges  are  designated  for  timing  and  ranging  (TR)  along  with 
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a  group  for  control  signaling  (C).  The  data  rate  is  2.048  Mbits/sec  and  the  entire 
payload  signal  is  confined  to  a  bandwidth  of  approximately  425  kHz  (see  markers 
on  the  spectral  plot).  The  bit  error  rate  (BER)  was  measured  as  a  function  of  the 
signal-to-noise  ratio  (SNR)  for  both  16-QAM  and  64-QAM  constellations  (Figure  37). 
The  measured  results  closely  track  the  calculated  (solid  line)  values.  The  calculations 
are  based  on  standard  QAM  performance  modeling  in  the  presence  of  Gaussian  noise 
[78]  and  the  SNR  is  measured  digitally  at  the  decision  statistic  after  equalization. 

Figure  38  shows  the  spectrum  of  the  signal  (which  appears  white  over  the  se¬ 
lected  frequency  band)  and  inserted  radio  frequency  interference  (RFI).  The  SNR  for 
each  tone  was  measured  over  the  band  for  different  power  levels  of  RFI  {Pingress  = 
—47,  —35,  —22  dBm)  with  the  signal  power  per  tone  held  constant  (Ptone  =  —41  dBm) 
Figure  39  graphs  the  measured  results,  demonstrating  how  the  DWMT  waveform  pro¬ 
vides  excellent  isolation  of  the  ingress  noise.  Even  under  the  worst  case  RFI  level, 
only  5  tones  (approximately  20  kHz  of  the  spectrum)  were  degraded  by  the  noise. 
The  remaining  tones  outside  of  this  band  attain  the  same  SNR  whether  or  not  the 
RFI  is  present. 

After  injecting  the  RFI  directly  in-band  as  shown  in  Figure  38,  a  group  of  tones 
were  reallocated  to  a  clean  portion  of  the  channel  and  measurements  were  made  as  a 
function  of  the  frequency  difference  from  the  center  of  the  closest  DWMT  tone  to  the 
center  of  the  RFI.  The  spectrum  of  the  corresponding  signal  with  64  QAM  modulation 
levels  and  ingress  noise  (RFI)  spectra  are  depicted  in  Figure  40.  With  no  RFI,  the 
SNR  for  each  DWMT  tone  is  about  35  dB.  Note  that  the  BER  was  measured  with  the 
system  subject  to  the  noise  levels  in  Figure  40  (signal/ingress  frequency  difference  is 
14  kHz).  No  errors  were  recorded  over  a  24  hour  period,  indicating  a  BER  less  than 

1  X  10~^^.  As  the  ingress  noise  was  moved  closer  to  the  signal,  the  BER  degrades  as 
shown  in  Figure  41.  The  measured  results  again  clearly  demonstrate  the  ability  of 
the  DWMT  system  to  isolate  the  ingress  noise.  Changing  the  ingress  frequency  by 

2  kHz  reduces  the  BER  from  10“^  to  10“^.  As  the  power  level  of  the  ingress  noise 
is  increased,  the  corresponding  frequency  separation  between  the  DWMT  signal  and 
the  ingress  must  also  increase.  This  is  a  direct  result  of  the  ingress  spectrum  that 
consists  of  a  main  lobe  and  a  series  of  side  lobes. 

In  conclusion,  wavelet  filter  banks  have  been  shown  both  theoretically  and  exper¬ 
imentally  to  add  a  new  dimension  to  multicarrier  modulation.  The  use  of  overlap  (as 
measured  by  the  genus  parameter  g)  provides  a  means  for  tradeoff  between  subchannel 
isolation  and  system  latency  and  computational  requirements.  The  subchannel  isola¬ 
tion  provided  by  wavelet  modulation  yields  increased  robustness  and  fiexibility  with 
regard  to  narrowband  noise.  Fast  algorithms  for  wavelet  computation  map  directly  to 
the  multicarrier  modulation  application,  minimizing  the  additional  cost  due  to  the  use 
of  an  overlapped  transform.  DWMT  is  a  very  promising  technique  for  high-bandwidth 
transmission  over  both  twisted-pair  copper  wire  and  hybrid  fiber-coax  wirelines  in  the 
“last  mile”  of  network  connections  to  home,  school,  and  business. 
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Figure  34;  Schematic  diagram  of  a  Discrete  Wavelet  Multitone  (DWMT)  system. 


Figure  35:  Experimental  baseband  DWMT  setup  with  1  head  end  modem  and  3 
remote  modems. 
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Figure  36:  Measured  spectrum  of  a  downstream  DWMT  signal  using  64  QAM  mod¬ 
ulation  for  the  payload  tones.  In  addition  to  the  tones  designated  for  payload  (P), 
additional  tones  are  used  for  ranging  (R)  and  control  (C). 


Figure  37:  Measured  bit  error  rate  (BER)  vs.  SNR  for  16-QAM  and  64-QAM  trans¬ 
mission  with  DWMT.  Solid  lines  are  the  theoretical  predictions. 
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Figure  38:  Spectral  plot  of  a  DWMT  signal  (appears  as  a  white  spectrum)  with  RFI 
present. 


Figure  39:  Measured  SNR  for  each  tone  with  RFI  at  power  levels  of  -47  dBm,  -35 
dBm,  and  -22dBm.  Even  with  very  high  RFI  levels,  only  5  tones  are  degraded  by  the 


noise. 
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Figure  40:  Spectral  plot  of  a  DWMT  signal  with  tones  reallocated  to  avoid  the  RFI. 
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Figure  41:  Measured  BER  for  the  noise  conditions  shown  in  Figure  40  for  three 
different  levels  of  RFI  and  64  QAM  signalling.  As  center  frequency  of  the  RFI  is 
varied  to  bring  the  RFI  within  the  band  of  a  tone,  the  BER  degrades. 


6.2  Image  Compression 

Image  data  compression  has  been  one  of  the  first  and  foremost  applications  of  wavelet 
methods.  We  continued  Aware’s  efforts  in  wavelet-based  image  compression  as  part 
of  the  higher  rank  wavelet  project.  Our  work  involved  applying  the  new  rank  M 
wavelet  transforms  to  image  compression,  fine-tuning  the  quantization  algorithms, 
and  applying  wavelet  techniques  to  the  compression  of  new  types  of  data,  including 
seismic  data  and  multispectral  LANDSAT  data. 

6.2.1  Review  of  wavelet-based  image  compression 

We  begin  with  a  review  of  wavelet-based  image  compression.  A  typical  lossy  com¬ 
pression  algorithm  consists  of  three  basic  steps:  a  reversible  lossless  transform,  a 
lossy  quantizer  and  a  standard  lossless  encoder.  The  transform  step  may  consist  of  a 
Discrete  Cosine  Transform  (DCT)  as  in  the  ANSI  standard  JPEG  compression  algo¬ 
rithm  [75],  or  a  wavelet  transform  as  in  Aware’s  AccuPress  production  software,  and 
in  the  research  reported  here.  The  energy  compaction  properties  of  the  transform 
domain  representation  is  crucial  to  the  overall  success  of  a  compression  algorithm. 
While  the  DCT  provides  a  compact  data  representation  composed  of  stationary  har¬ 
monics,  wavelet  transforms  provide  a  superior  representation  for  data  composed  of 
sharp  edges  and  other  transient  characteristics  as  well  as  broad  smooth  areas.  The 
quantization  step,  which  accounts  for  all  of  the  “loss”  in  the  compression,  seeks  to 
represent  the  more  significant  components  of  the  transform  representation  with  pro¬ 
portionally  more  accuracy  or  bits  than  the  less  significant  components.  The  ability  of 
the  wavelet  transform  to  compactly  segregate  the  important  parts  of  seismic  signals 
allows  the  quantizer,  for  a  given  bitrate,  to  cause  significantly  less  “loss”  than  would 
be  the  case  for  other  representations  such  as  a  Fourier  transform  [64].  Because  the 
quantizer  output  is  put  through  a  lossless  entropy  coder  to  create  the  final  compressed 
bitstream,  we  employ  an  entropy-constrained  quantization  technique.  The  quantized 
data  is  then  passed  through  a  combined  Huffman  and  zero-run-length  encoding  to 
approach  the  true  entropy  of  the  quantized  bitstream. 

6.2.2  New  wavelet  transforms 

A  significant  piece  of  the  work  carried  out  in  this  project  was  the  development  of 
new  families  of  wavelet  transforms  (the  rank  M  wavelets).  We  applied  several  of 
the  new  wavelet  constructions  to  image  compression  using  Aware’s  standard  wavelet 
compression  algorithm;  results  were  briefly  quoted  in  sections  4.4.4  and  4.5.2.  Further 
details  can  be  found  in  references  [45],  [48],  and  [44]. 
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6.2.3  Improved  quantization  techniques 

Our  efforts  on  improving  (lossy)  quantization  were  twofold:  an  examination  of  so- 
called  “dead  zones”  for  wavelet  scalar  quantizers,  and  the  development  of  an  opera¬ 
tional  rate-distortion-optimal  uniform  quantization  scheme. 

Empirical  evidence  [33]  has  suggested  that  the  subjective  and  objective  performace 
of  a  uniform  quantizer  may  be  improved  in  transform  coding  applications  by  allowing 
the  center  decision  region  in  a  uniform  quantizer  (also  known  as  the  zero  bin)  to  be 
wider  than  the  other  decision  regions  by  a  constant  factor.  In  our  notation,  the  zero 
bin  width  is  given  by 

Ao  =  /?A 

where  A  is  the  uniform  quantizer  bin  width. 

The  expanded  zero  quantization  bin  is  sometimes  called  the  dead  zone.  It  is 
thought  that  the  perceptual  improvement  afforded  by  an  expanded  zero  bin  is  due  to 
the  fact  that  the  contribution  of  basis  functions  with  small  coefficients  is  underem¬ 
phasized  rather  than  exaggerated.  While  this  does  not  minimize  the  mean  squared 
error,  it  improves  perceptual  quality  because  small  coefficients  increased  by  amplitude 
quantization  manifest  themselves  perceptually  as  noise.  With  an  expanded  dead  zone 
quantizer,  more  of  the  extremely  small  coefficients  are  quantized  to  zero,  yielding  a 
more  visually  pleasing  reconstruction  and  aiding  subsequent  zero-run-length  coding. 

We  conducted  experiments  with  varying  dead  zone  widths  on  a  variety  of  image 
data,  ranging  from  natural  images  to  overhead  satellite  images  to  seismic  datasets. 
We  found  that  using  an  expanded  center  bin  yielded  an  increase  in  peak  signal-to- 
noise  ratio  (PSNR)  between  0.2  to  0.5  dB.  For  natural  and  satellite  imagery,  the 
maximally  effective  dead  zone  width  (around  /5  =  1.8)  yielded  a  PSNR  improvement 
of  about  0.3  dB.  Both  the  location  of  the  peak  and  the  amount  of  improvement  seem 
to  be  typical  for  natural  images  such  as  the  NITF  series  at  this  compression  ratio. 
As  the  compression  ratio  increase,  the  peak  PSNR  improvement  occurred  for  larger 
values  of  /?.  The  most  significant  difference  between  the  seismic  data  and  the  natural 
images  was  that  the  peak  SNR  improvement  seemed  to  occur  at  a  much  lower  /?  for 
the  seismic  data.  In  order  to  compress  both  types  of  image  data  effectively,  we  chose 
a  fixed  (3  of  1.4,  which  gave  reasonable  performance  improvement  for  both  data  types. 
Further  details  appear  in  [17]. 

In  addition  to  this  examination  of  dead-zone  quantizers,  we  developed  an  op¬ 
erational  rate-distortion-optimized  quantizer  based  on  the  methods  of  Shoham  and 
Gersho  [93],  applied  to  the  wavelet  transform  domain.  It  has  been  observed  [108],  [64] 
that  histograms  of  the  subband  outputs  of  the  wavelet  transform  applied  to  image 
data  obey  a  generalized  Gaussian  distribution.  Beginning  with  this  assumption,  we 
employed  an  entropy-constrained  uniform  quantizer  with  a  mid-rise  transfer  curve. 

The  bit- allocation  algorithms  employed  [93],  [80]  require  knowledge  of  the  rate- 
distortion  characteristics  of  the  quantized  data.  Since  our  quantizer  is  followed  by  an 
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entropy  encoder,  we  make  the  assumption  that  the  source  will  be  compressed  down 
to  its  entropy  level.  For  this  reason,  we  examine  the  entropy  of  various  quantized 
sources  versus  the  distortion  rather  than  rate  versus  distortion.  We  also  chose  uniform 
quantizers  as  opposed  to  the  pdf-optimized  Lloyd-Max  quantizers  ([108],  [55]). 

First,  we  compute  the  entropy  of  the  quantized  source  as  a  function  of  the  bin 
width  A  and  the  generalized  Gaussian  exponent  a.  The  entropy  is  defined  to  be 
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H  Pi  log2  Pi 
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that  is,  a  sum  of  integrals  of  the  probability  distribution  of  the  source  over  the  i 
decision  region.  Using  these  expressions,  we  can  easily  compute  the  value  of  the 
entropy  as  a  function  of  the  stepsize  A  and  the  subband  variance  cr. 

The  distortion  is  defined  to  be  the  mean  squared  error  introduced  by  the  quanti¬ 
zation  process.  It  is  clear  that  as  the  step  size  becomes  finer  and  finer,  the  approx¬ 
imation  is  increasingly  better,  thus  yielding  smaller  distortion.  On  the  other  hand, 
the  more  finely  quantized  the  data  is,  the  more  bits  are  needed  to  represent  it.  This 
is  the  fundamental  tradeoff  involved  in  the  lossy  transmission  of  information  the 
rate/distortion  tradeoff. 

The  distortion  due  to  quantization  in  a  given  subband  is 
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Since  the  quantization  bins  form  a  partition  of  the  real  line,  we  may  write  the  integral 
as  a  summation  of  integrals  over  each  quantization  bin  separately  plus  the  integral 
over  the  tails  of  the  distribution. 
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where  Vi  is  the  i-th  reconstruction  level.  Using  the  component  parts  of  this  equa¬ 
tion,  we  can  compute  the  overall  distortion  introduced  by  the  quantizer.  Having 
now  parameterized  both  the  entropy  H  and  the  distortion  D  by  the  binwidth,  sub¬ 
band  variance,  and  Gaussian  exponent,  we  can  now  compute  parameterized  tables  of 
entropy  versus  distortion  for  use  in  quantization  and  bit  allocation. 

Given  such  precomputed  rate-distortion  tables  for  a  fixed  set  of  generalized  Gaus¬ 
sian  exponents,  at  runtime  we  measure  the  subband  variance  and  kurtosis  to  de¬ 
termine  the  best-fit  Gaussian  exponent  to  each  subband.  This  determines  which 
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rate-distortion  tables  are  used  to  determine  the  optimal  bit  allocation  among  the 
subbands,  using  the  operational  rate-distortion  methodology  of  [93],  which  follows 
the  convex  hull  of  the  rate-distortion  curve. 

The  resulting  quantizer  yields  significantly  superior  performance  for  wavelet  com¬ 
pression,  when  compared  with  simple  scalar  quantizers  such  as  that  of  the  FBI’s 
Wavelet  Scalar  Quantization  specification  [7].  Further  details  are  given  in  [18]. 

6.2.4  Seismic  data  compression 

In  addition  to  improvement  of  our  basic  wavelet  transform  compression  algorithm,  we 
applied  wavelet  compression  to  several  new  types  of  data.  One  of  the  most  promising 
developments  has  been  the  application  to  seismic  data  compression.  In  the  seismic 
regime  we  tested  a  variety  of  wavelet  bases  and  tree  structures  on  a  variety  of  seismic 
data  types.  After  a  brief  review  of  the  specific  areas  of  experimentation,  we  explore 
our  prospects  and  successes  in  commercialization  of  wavelet-based  seismic  data  com¬ 
pression. 

One  of  our  first  areas  of  experimentation  was  to  try  different  rank  2  basis  func¬ 
tions  (filters)  and  different  multiresolution  trees  on  sample  seismic  datasets.  The 
Daubechies  biorthogonal  7/9-tap  filters  [11]  were  compared  to  Daubechies  6,8,10,  and 
12  tap  orthogonal  filters.  The  biorthogonal  7-9  tap  filters  provided  up  to  a  several  dB 
increase  in  pSNR  (peak  signal-to-noise  ratio)  over  any  of  the  Daubechies  orthogonal 
filters  at  compression  ratios  of  between  20:1  and  80:1.  This  is  consistent  with  the 
results  of  Macq  and  Mertes  on  natural  images  [63]. 

We  also  explored  the  use  of  new  families  of  rank  M  wavelets;  in  particular,  the 
rank  M  linear-phase  orthogonal  wavelets  discussed  in  section  4.5.2.  Compression 
results  (based  on  the  size  of  the  entropy-coded  bitstream)  are  shown  in  Table  3.  We 
compared  a  5-level  Mallat  tree  based  on  the  rank  2  Daubechies  (7,9)-tap  biorthogonal 
filter  pair  [11]  with  a  3-level  Mallat  tree  [64]  based  on  the  rank  4  12-tap  GenLOT 
shown  in  Figure  25.  The  rank  4  wavelet  offered  superior  or  comparable  performance 
across  all  compression  ratios;  peak  SNR’s  were  approximately  1  dB  greater  for  the 
rank  4  linear- phase  wavelet. 


8:1 

16:1 

32:1 

pSNR 

Max 

pSNR 

Max 

pSNR 

Max 

M  =  2 

51.6 

1101 

40.7 

4262 

32.7 

12544 

M  =  4 

52.6 

1025 

41.3 

4298 

33.7 

12260 

Table  3:  Peak  SNR  and  maximum  errors  for  compression  of  2-d  seismic  data  example 
(Mallat  tree,  M  —  2  (7,9)-tap  pair  and  M  =  4  GenLOT). 

We  also  tested  a  number  of  non-Mallat  transform  trees  to  evaluate  the  merit  of  us- 
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ing  hardcoded  trees  for  a  seismic-specific  wavelet  compression  algorithm.  The  recently 
released  Federal  Bureau  of  Investigation  specification  for  fingerprint  compression  has 
such  a  hardcoded  tree  which  was  found  to  be  optimal  for  fingerprints  (Figure  27). 
We  evaluated  several  non-Mallat  trees  on  four  very  different  types  of  seismic  data. 
Our  initial  conclusion  is  that  the  wide  variation  in  seismic  data  types  makes  the  use 
of  a  single  fixed  non-Mallat  tree  impractical.  However,  a  fast  adaptive  tree  could  be 
successful  for  seismic  applications. 

The  efficiency  of  a  standard  two-dimensional  wavelet  transform  for  data  which  is 
moderately  rectangular  was  found  to  be  generally  adequate.  Seismic  data  tends  to 
be  grouped  into  logical  frames  with  approximate  dimensions  of  240  columns  by  2000 
rows.  We  compared  a  standard  Mallat  tree  with  an  inital  rowwise  transform  followed 
by  a  Mallat  tree  on  the  low  pass  output  and  found  the  results  to  be  comparable. 

Seismic  data  can  be  grouped  into  two  broad  categories,  stacked  and  pre-stack. 
Stacked  data  contains  a  large  percentage  of  horizontally  aligned  energy  caused  by 
the  generally  fiat  structure  of  the  earth’s  sedimentary  rock  layers.  Pre-stack  data  is 
composed  largely  of  hyperbolically  aligned  energy  due  to  the  geometry  of  the  seismic 
data  acquisition  process.  In  theory  the  horizontally  aligned  arrivals  should  be  easier 
to  compress,  but  in  practice  the  hyperbolically  aligned  seismic  arrivals  were  preserved 
as  well  as  the  horizontally  aligned  arrivals  after  passing  through  the  same  wavelet- 
transform  based  compression  algorithm.  This  was  the  opinion  of  the  seismic  experts 
who  examined  the  data. 

Another  area  of  algorithm  development  lay  in  the  use  of  companding.  Condition¬ 
ing  or  companding  the  histogram  of  data  prior  to  compression  is  a  known  method  for 
improving  compression  performance  [55].  This  is  due  to  improved  quantizer  perfor¬ 
mance  because  the  companding  operation  has  in  some  way  caused  the  quantizer  to 
preferentially  preserve  amplitudes  which  are  important  to  the  application  at  hand. 
The  largest  amplitudes  in  a  seismic  data  set  are  usually  the  least  important  compo¬ 
nents  of  the  data  and  therefore  should  not  be  preserved  at  the  expense  of  other  more 
important  components.  While  it  is  not  possible  to  simply  compress  the  reciprocal 
of  the  seismic  data,  it  does  help  to  condition  the  data  prior  to  compression  to  help 
boost  small  but  important  features.  We  experimented  with  standard  logarithm  based 
companding  curves  and  found  moderate  improvement  at  low  compression  ratios  but 
unacceptable  results  for  compression  ratio  greater  than  15:1.  Much  better  results 
have  been  achieved  by  allowing  the  seismic  data  processing  specialist  to  do  this  con¬ 
ditioning  based  on  some  physical  model  of  the  earth  prior  to  compression.  Specific 
characteristics  of  the  exploration  regime  can  be  used  in  the  conditioning  of  the  data, 
rather  than  a  simplistic  blind  companding  curve.  Initial  results  indicate  that  con¬ 
ditioning  of  seismic  data  is  an  important  step  prior  to  compression,  particularly  for 
pre-stack  data. 

We  see  considerable  commercial  potential  for  wavelet-based  seismic  data  com¬ 
pression.  A  typical  seismic  acquisition  boat  collects  between  one  and  ten  gigabytes 
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of  seismic  data  per  day.  Historically,  the  seismic  industry  has  been  the  initial  con¬ 
sumer  very  large  mass  storage  and  data  transmission  products  and  continues  to  do 
so.  Aware  has  worked  with  a  number  of  companies  involved  in  actual  oil  and  gas  ex¬ 
ploration  as  well  as  in  data  acquisition  and  data  processing.  Lossy  data  compression 
is  a  novel  concept  to  the  entire  seismic  industry  and  has  been  quite  well  received.  We 
have  had  experts  in  seismic  data  processing  state  that  the  compression  results  at  10:1 
are  virtually  ’’lossless”  and  20:1  is  sufficient  for  the  vast  majority  of  processing.  We 
have  successfully  integrated  compression  into  data  transfer  systems  which  use  the  In¬ 
marsat  communications  satellites.  These  systems  allow  technicians  on  seismic  boats 
to  send  samples  of  seismic  data  via  Inmarsat  back  to  a  home  office  where  experts 
can  examine  the  data  and  determine  whether  changes  need  to  be  made  to  the  boat’s 
acquisition  plans.  While  the  file  transfer  systems  previously  existed,  it  was  virtually 
impossible  to  send  seismic  data  due  to  the  cost  of  the  satellite  link.  Data  compression 
has  had  an  enormous  impact  on  costs,  enabling  the  practical  economic  transmission  of 
seismic  data  from  acquisition  vessels  in  near  real  time.  We  are  also  working  on  using 
compression  to  reduce  seismic  mass  storage  costs  in  on-shore  processing  centers.  Our 
close  contact  with  the  seismic  industry  has  allowed  the  algorthimic  efforts  detailed 
above  to  be  reviewed  by  a  demanding  and  knowledgeable  audience.  Further  details 
appear  in  the  publications  [14],  [83]. 

6.2.5  Compression  of  Multispectral  (LANDSAT)  Data 

As  part  of  the  contract  work  effort,  we  also  extended  lossy  wavelet  image  compression 
methods  to  12-band  multispectral  LANDSAT  satellite  imagery.  Multispectral  and 
hyperspectral  images,  usually  gathered  by  satellite,  consist  of  many  two-dimensional 
images.  Each  two-dimensional  image  corresponds  to  a  narrow  spectral  band.  The 
LANDSAT  data  we  worked  with  had  12  such  spectral  planes  or  bands;  next-generation 
satellites  capture  hyperspectral  images  with  60  to  100  planes.  Each  band  is  usually 
about  1000  X  1000  pixels,  with  8  bits  per  pixel.  There  is  significant  correlation  among 
the  spectral  bands;  exploiting  this  correlation  is  key  to  successful  multispectral  data 
compression. 

In  practice,  the  multispectral  data  will  be  gathered  and  compressed  on  board  a 
satellite,  then  transmitted  to  a  ground  station  for  storage,  decompression,  and  re¬ 
trieval.  The  satellite  will  have  strict  limitations  on  power  consumption  and  hence  the 
possible  complexity /memory  requirements  of  the  compression  algorithm,  and  must 
operate  in  real  time.  In  contrast,  the  ground  station  will  have  practically  nonexistent 
power  limitations,  and  the  received  or  stored  data  will  often  be  browsed  in  a  non-real¬ 
time  setting.  In  view  of  this  asymmetrical  system,  several  requirements  appear.  The 
compression  algorithm  must  not  have  too  high  a  complexity.  It  must  also  be  of  high 
quality,  for  the  satellite  cannot  go  back  and  recapture  data  that  has  been  found  to  be 
of  low  fidelity.  Furthermore,  a  scheme  which  lends  itself  to  searching  and  browsing  is 
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preferable. 

A  wavelet  approach  wins  on  several  of  these  counts.  Wavelet  compression  has 
consistently  proven  itself  to  be  of  higher  quality  than  comparable  DCT-based  algo¬ 
rithms  [16],  while  of  only  slightly  greater  complexity.  Because  the  wavelet  transform 
is  a  multiresolution  representation,  it  offers  several  advantages  for  LANDSAT  appli¬ 
cations.  A  few  coarse-scale  coefl&cients  can  efficiently  describe  broad  subregions  of  a 
large  satellite  image  while  yielding  high  compression  ratios.  This  is  in  sharp  distinc¬ 
tion  to  the  blocking  artifacts  and  limited  effective  compression  ratios  of  block-DCT 
algorithms  such  as  the  JPEG  standard  [75].  Wavelet  data  structures  also  naturally 
lend  themselves  to  multiresolution  image  browsing,  allowing  the  user  to  call  up  a 
coarse-scale  representation  of  an  image  and  only  request  finer  detail  when  he  actually 
needs  it. 

Our  experiments  involved  the  exploration  of  several  alternatives  in  the  design  of 
a  compression  algorithm  for  multispectral  data,  within  the  constraint  of  employing  a 
wavelet  transform  in  each  of  the  two-dimensional  bands  of  the  multispectral  dataset. 
The  LANDSAT  test  imagery  was  of  varying  sizes,  roughly  900  x  1100  pixels  of  12 
bands  by  8  bits.  We  applied  a  wavelet  transform  in  the  spatial  dimension,  consisting 
of  a  six  level  Mallat  tree  of  Daubechies’  rank  2  linear-phase  (7,9)-tap  biorthogonal 
wavelets  [11],  yielding  a  nonexpansive  transform  of  the  data.  We  then  considered  the 
following  alternatives: 

•  Independent  wavelet  compression  processing  of  each  spectral  band. 

•  Independent  wavelet  transforms  of  each  band,  followed  by  global  bit  allocation 
(quantizing  all  subbands  of  all  spectral  bands  at  once). 

•  The  use  of  a  Discrete  Cosine  Transform  (DCT)  in  the  12-band  spectral  dimen¬ 
sion  as  well  as  wavelet  transforms  spatially,  followed  by  global  bit  allocation. 

•  The  use  of  the  Karhunen-Loeve  Transform  (KLT)  in  the  12-band  spectral^  di¬ 
mension  as  well  as  wavelet  transforms  spatially,  followed  by  global  bit  allocation. 

The  first  two  alternatives  were  considered  as  a  baseline,  against  which  to  measure 
the  elTectiveness  of  the  spectral  DCT  and  KLT  transforms  in  exploiting  cross-band 
redundancy.  The  only  difference  between  the  first  two  approaches  is  in  the  bit  alloca¬ 
tion,  where  the  first  employed  distinct  bit  allocations  for  each  band,  while  the  second 
approach  allocated  bits  across  all  the  subbands.  Not  surprisingly,  the  two  approaches 
displayedd  very  similar  performance. 

We  then  considered  the  latter  pair  of  alternatives,  performing  either  a  Discrete 
Cosine  Transform  or  a  Karhunen-Loeve  Transform  across  the  12  spectral  bands.  The 
KLT,  which  is  the  transform  that  diagonalizes  the  cross-correlation  matrix  of  the 
data,  is  precisely  that  block  transform  that  yields  optimal  energy  compaction  of  the 
data.  Because  of  this  data-dependent  optimality,  the  KLT  will  do  the  best  possible 
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job  of  transforming  the  cross-band  data  for  subsequent  quantization  and  compres¬ 
sion.  However,  we  considered  the  DOT  as  well  because  system  requirements  may 
dictate  the  use  of  a  fixed  (non-data-dependent)  transform  with  fast  algorithms  for 
on-board  compression  processing.  The  DCT  is  widely  used,  because  it  both  posesses 
fast  computational  algorithms  and  serves  as  an  approximation  to  the  KLT  for  AR(1) 
processes  [82]. 

Rate-distortion  curves  for  the  various  approaches  across  a  range  of  bitrates  from 
.25  bits/pixel  to  2  bits/pixel  (compression  ratios  ranging  from  32:1  to  4:1)  are  shown 
in  Figure  42,  while  compressed/decompressed  images  are  shown  as  Figure  43.  Observe 
the  significant  decrease  in  distortion  when  one  moves  from  compression  based  on  a 
spatial-only  transform  to  the  addition  of  a  cross-band  transform  such  as  the  DCT  or 
KLT.  Interestingly,  the  DCT  yielded  75%  of  the  performance  gain  achieved  by  going 
from  no  spectral  transform  to  the  full-blown  KLT. 


Figure  42:  Rate-distortion  curves  for  the  four  multispectral  compression  algorithms. 

In  conclusion,  these  experiments  (reported  in  [41])  confirm  the  utility  of  wavelet 
transforms  for  multispectral  image  compression,  and  the  importance  of  exploiting 
cross-band  correlation.  We  also  validated  the  use  of  a  DCT  as  a  stand-in  for  the  KLT 
in  exploiting  cross-band  correlations  in  a  computationally  efficient  manner.  Areas  for 
future  work  include  trying  simple  (e.g.  Haar)  wavelet  transforms  across  the  bands, 
employing  zerotree  structures  [92]  in  the  quantization  and  coding,  and  trying  out  ideas 
from  video  compression  like  motion  estimation  and  compensation  for  the  cross-band 
processing. 
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Figure  43:  See  attached  page. 


6.3  Sonar  Data  Compression  in  Real  Time 


During  1993  Aware  took  part  in  an  AntiSubmarine  Warfare  simulation/exercise  using 
the  Internet  that  was  sponsored  by  ARPA’s  Maritime  Systems  and  Technology  Office 
[25],  Our  role  was  to  provide  high  performance  data  compression  to  reduce  inter-site 
communications  bandwidth  requirements  in  a  distributed  simulation  environment. 
Aware  has  already  developed  robust  CD-quality  real-time  audio  compression  software 
[98,  99]  which  is  based  on  rank  M  wavelets  and  multirate  filtering.  A  critical  part  of 
sonar  compression  project  was  to  tailor  the  basis  choice  to  sonar  signals,  rather  than 
CD  audio.  The  WaveTool  software  proved  invaluable  in  this  task. 

Given  sample  sonar  data,  the  aim  was  to  find  a  wavelet/ subband  decomposition 
which  provided  maximum  energy  compaction,  i.e.  localized  the  signal  energy  in  as 
few  subbands  as  possible.  We  used  a  well-known  bit  allocation  formula  that  optimizes 
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to  achieve  an  overall  bitrate  b.  We  decided  to  restrict  ourselves  to  the  use  of  cosine- 
modulated  perfect  reconstruction  filter  banks,  because  of  the  available  DCT-based 
fast  algorithms  for  their  computation,  as  well  as  the  superiority  of  off-the-shelf  filter 
design  methods.  We  ran  several  sonar  test  signals  through  the  Wavetool,  testing  filter 
banks  of  rank  4,  8,  and  16.  In  each  case,  we  analyzed  the  input  signal,  wrote  the 
subband  data  out  to  file,  quantized  using  the  scheme  (43),  and  then  read  the  quantized 
data  back  in  to  the  Wavetool  and  synthesized  a  reconstructed  signal.  Using  the  error 
measurement  features  of  the  Wavetool,  we  were  able  to  evaluate  mean-squared  error 
and  maximum  error,  and  found  the  rank  8  filter  banks  to  be  optimal  for  the  signals 
in  question.  We  then  varied  the  overlap  of  the  M  =  8  filter  bank  and  decided  on 
an  overlap  =  4  design.  Longer  overlaps  led  to  diminishing  returns  at  a  nonzero 
computational  cost.  The  sonar  data  to  be  compressed  was  classified  as  active  and 
passive.  The  passive  data  was  principally  lowpass,  and  using  the  Wavetool,  we  found 
that  we  could  attain  further  compaction  of  the  passive  signals  by  iterating  on  the 
lowest  channel  of  the  M  8  split  with  a  rank  2  split.  For  this  filter  we  used  a  64-tap 
high  stopband  attenuation  perfect-reconstruction  bank.  The  facilities  provided  by 
the  WaveTool  led  to  a  significant  reduction  in  algorithm  design  time  for  the  sonar 
data  compression  demonstration. 


6.4  Multiwavelet  Signal  Processing 

A  very  recent  development  in  wavelet  theory  is  the  class  of  multiwavelets  -  wave¬ 
let  decompositions  based  on  matrix  dilation  equations.  While  these  structures  have 
been  developed  by  mathematicians,  significant  gaps  in  application  of  multiwavelets 
presented  themselves  to  us.  We  have  addressed  the  computational  issues  of  pre¬ 
filtering  and  symmetric  extension  for  multiwavelet  signal  and  image  processing,  and 
applied  multiwavelet  filters  to  image  compression  and  denoising  [49],  [105],  [106].  We 
summarize  these  results  here. 

Multiwavelets  are  also  based  on  a  multiresolution  analysis,  however  one  that  is 
based  on  several  scaling  functions.  A  basis  for  the  coarse  approxmiation  space  Vq  is 
generated  by  translates  of  N  scaling  functions  4>i{t  —  k),  —  k),...,  4>N{i  —  k).  The 

vector  $(t)  =  [</>i(t),  . .  • ,  <^Ar(t)]^,  will  satisfy  a  matrix  dilation  equation  (analogous 
to  the  scalar  case) 

Ht)  =  '£C[k]^2t-k).  (44) 

k 

The  coefficients  C[k]  are  N  hy  N  matrices  instead  of  scalars. 

Associated  with  these  scaling  functions  are  N  wavelets  Wi(t), . .  ■ ,  satisfying 

the  matrix  wavelet  equation 

W(t)  =  ^D[A;]#(2t-fc).  (45) 

k 

Again,  W{t)  =  [r£;i(t),  . . . ,  WN{t)]'^  is  a  vector  and  the  D[k]  are  A"  by  iV  matrices. 

As  in  the  scalar  case,  one  can  find  the  conditions  of  orthogonality  and  approxi¬ 
mation  for  multiwavelets  [103,  104,  37,  76]. 


Figure  44:  Geronimo-Hardin-Massopust  pair  of  scaling  functions. 

A  very  important  multiwavelet  system  was  constructed  by  J.  Geronimo,  D.  Hardin, 
and  P.  Massopust  [32]  (see  [9]  for  another  early  multiwavelet  construction).  Their 
system  contains  the  two  scaling  functions  <^i(t),  </>2(t)  shown  in  Figure  44  and  the  two 
wavelets  Wi(^t),W2{t)  shown  in  Figure  45.  The  dilation  and  wavelet  equations  for  this 
system  have  four  coefficients: 

Hi)  =  [  =  C[0M2t)  +  C[l]H2t  -  1)  -b  C[2]$(2t  -  2)  +  C'[3]$(2t  -  3) , 

[  92{t) 
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Figure  45:  Geronimo-Hardin-Massopust  multi  wavelets. 


W{t):= 


Wi{t) 

W2it) 


D[0]$(2t)  +  D[im2t  -  1)  +  D[2M2t  -  2)  +  D[3]^2t  -  3) , 


Actual  values  for  the  coefficient  matrices  C\i]  and  D\i]  can  be  found  in  [105].  There 
are  four  remarkable  properties  of  the  Geronimo-Hardin-Massopust  scaling  functions. 

•  They  each  have  short  support  (the  intervals  [0,  1]  and  [0,  2]). 

•  Both  scaling  functions  are  symmetric,  and  the  wavelets  form  a  symmetric/anti¬ 
symmetric  pair. 

•  All  integer  translates  of  the  scaling  functions  are  orthogonal. 

•  The  system  has  second  order  of  approximation  (locally  constant  and  locally 
linear  functions  are  in  ho). 

Let  us  stress  that  a  scalar  system  with  one  scaling  function  cannot  combine  sym¬ 
metry,  orthogonality,  and  second  order  approximation.  Moreover,  a  solution  of  a 
scalar  dilation  equation  with  four  coefficients  is  supported  on  the  interval  [0,  3]! 
Other  useful  constructions  of  multiwavelets  are  given  in  [9,  30,  107,  112,  77,  61]. 

6.4.1  Multiwavelets  and  multirate  filter  banks 

Corresponding  to  each  multiwavelet  system  is  a  matrix- valued  multirate  filter  bank  or 
multifilter.  A  multiwavelet  filter  bank  [103]  has  “taps”  that  are  N  x  N  matrices  (we 
work  primarily  with  =  2).  The  principal  example  is  the  4-coefficient  symmetric 
multiwavelet  filter  bank  whose  lowpass  filter  was  reported  in  [32].  This  filter  is  given 
by  the  four  2x2  matrices  (7[A:].  The  corresponding  2-channel,  2x2  matrix  filter 
bank  operates  on  two  input  data  streams,  filtering  them  into  four  output  streams, 
each  of  which  is  downsampled  by  a  factor  of  2,  This  is  shown  in  Figure  46.  Each  row 
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of  the  multifilter  is  a  combination  of  two  ordinary  filters,  one  operating  on  the  first 
data  stream  and  the  other  operating  on  the  second.  For  example,  the  first  lowpass 
multiwavelet  filter  (7(0]  operates  as  co,o[^]  on  the  first  input  stream  and  Co,i[/s]  on  the 
second.  It  is  a  combination  of  the  Haar  filter  {1,1}  on  the  first  stream  and  the  unit 
impulse  response  on  the  second  stream. 


Figure  46:  A  multiwavelet  filter  bank,  iterated  once. 


The  matrix  filter  coefficients  must  satisfy  the  orthogonality  (“block-paraunitarity”) 
condition 

JV-l 

X]  C[k]  C[k  -  2lf  =  2So,i  I . 

A:=:0 


(46) 


The  lowpass  filter  C  and  highpass  filter  D  consist  of  coefficients  corresponding  to 
the  dilation  equation  (44)  and  wavelet  equation  (45).  But  in  the  multiwavelet  setting 
these  coefficients  are  n  by  n  matrices,  and  during  the  convolution  step  they  must  mul¬ 
tiply  vectors  (instead  of  scalars).  This  means  that  multifilter  banks  need  n  input  rows. 
We  developed  several  ways  to  produce  those  rows.  The  first  method,  an  oversampled 
scheme,  consisted  of  simply  repeating  each  row  of  data  twice  to  produce  the  two  in¬ 
put  rows  required.  This  method  yielded  good  results  at  denoising  one-dimensional 
signals,  but  performed  poorly  for  compression,  not  a  surprising  outcome.  Our  second 
method,  which  was  successful  at  both  compression  and  denoising,  exploited  the  ap¬ 
proximation  properties  of  the  multiwavelet  scaling  functions  to  prefilter  each  row  of 
input  data,  resulting  in  two  new  rows  with  half  the  number  of  data  points  for  input  to 
the  multifilter.  This  later  method  had  the  advantage  of  preserving  critical  sampling, 
so  that  an  input  dataset  with  N  points  remained  of  size  N  after  prefiltering.  Another 
advantage  of  this  approximation-based  preprocessing  method  is  that  it  fits  naturally 
with  symmetric  extension  for  multiwavelets  (discussed  below).  In  other  words,  if  we 
symmetrically  extend  a  finite  length  signal  f[n]  at  its  boundaries  and  implement  the 
approximation  formulas,  then  the  two  rows  output  by  the  preprocessor  will  have  the 
appropriate  symmetry. 

In  the  setting  of  purely  two-dimensional  signal  processing,  we  described  an  addi¬ 
tional  algorithm  for  multi  wavelet  filtering  (two  rows  at  a  time),  and  developed  a  new 
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family  of  multiwavelets  (the  constrained  pairs)  that  is  well-suited  to  this  two-row-at- 
a-time  filtering.  Further  details  appear  in  [105]. 

In  practice  all  signals  have  finite  length,  so  we  must  devise  techniques  for  filtering 
such  signals  at  their  boundaries.  There  are  two  common  methods  for  filtering  at  the 
boundary  that  preserve  critical  sampling.  The  first  is  circular  periodization  (periodic 
wrap)  of  the  data.  This  method  introduces  discontinuities  at  the  boundaries;  how¬ 
ever,  it  can  be  used  with  almost  any  filter  bank.  The  second  approach  is  symmetric 
extension  of  the  data.  Symmetric  extension  preserves  signal  continuity,  but  can  be 
implemented  only  with  linear-phase  (symmetric  and/or  antisymmetric)  filter  banks 
[95,  12,  56,  15].  We  have  developed  symmetric  extension  for  linear-phase  multiwavelet 
filters,  such  as  the  Geronimo-Hardin-Massopust  multifilters.  This  has  proven  useful 
for  image  compression  applications. 

Recall  the  basic  problem:  given  an  input  signal  f[n]  with  N  samples  and  a  linear- 
phase  (symmetric  or  antisymmetric)  filter,  how  can  we  symmetrically  extend  /  before 
filtering  and  downsampling  in  a  way  that  preserves  the  critically  sampled  nature  of 
the  system?  The  possibilities  for  such  an  extension  have  been  enumerated  in  [15]. 
Depending  on  the  parity  of  the  input  signal  (even-  or  odd-length)  and  the  parity 
and  symmetry  of  the  filter,  there  is  a  specific  non-expansive  symmetric  extension  of 
both  the  input  signal  and  the  subband  outputs.  For  example,  an  even-length  input 
signal  passed  through  an  even-length  symmetric  lowpass  filter  should  be  extended  by 
repeating  the  first  and  last  samples,  i.e.,  a  half-sample  symmetric  signal  is  matched 
to  a  half-sample-symmetric  filter.  Similarly,  when  the  lowpass  filter  is  of  odd  length 
(whole-sample-symmetry),  the  input  signal  should  be  extended  without  repeating  the 
first  or  last  samples. 

Each  row  of  the  GHM  multifilter  is  a  linear  combination  of  two  filters,  one  for 
each  input  stream.  One  filter  (applied  to  the  first  stream)  is  of  even  length;  the  sec¬ 
ond  is  of  odd  length.  Thus  we  extend  the  first  stream  using  half-sample-symmetry 
(repeating  the  first  and  last  samples)  and  extend  the  second  stream  using  whole- 
sample-symmetry  (noi  repeating  samples).  Then,  when  synthesizing  the  input  signal 
from  the  subband  outputs,  we  must  symmetrize  the  subband  data  differently  depend¬ 
ing  on  whether  it  is  going  into  an  even-  or  odd-length  filter.  This  approach  has  been 
extended  to  multifilters  lacking  linear-phase  symmetry  as  well.  Details  appear  in 
[105];  the  upshot  is  that  we  have  obtained  a  non-expansive  transform  of  finite-length 
input  data  which  behaves  well  at  the  boundaries  under  lossy  quantization. 

6.4.2  Denoising  by  soft  thresholding 

We  have  compared  the  numerical  performance  of  GHM  and  constrained  multiwavelets 
with  Daubechies  D4  scalar  wavelets.  D4  wavelets  were  chosen  because  they  have  two 
vanishing  moments,  are  orthogonal  and  have  four  coefficients  in  the  dilation  equation 
—  exactly  like  the  GHM  and  constrained  pairs.  We  perform  these  comparisons  in 
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two  standard  wavelet  applications:  signal  denoising  and  data  compression.  First  we 
discuss  the  denoising. 

Suppose  that  a  signal  of  interest  /  has  been  corrupted  by  noise,  so  that  we  observe 
a  signal  g: 


g[n]  =  f[n]  +  az[n], 

where  z[n\  is  unit- variance,  zero-mean  Gaussian  white  noise.  What  is  a  robust  method 
for  recovering  /  from  the  samples  g[rii\  as  best  as  possible?  Donoho  and  Johnstone 
[28,  29]  have  proposed  a  solution  via  wavelet  shrinkage  or  soft  thresholding  in  the 
wavelet  domain.  Wavelet  shrinkage  works  as  follows: 

1.  Apply  the  cascade  algorithm  to  get  the  wavelet  coefficients  corresponding 
to  g[nY 

2.  Choose  a  threshold  tn  =  y^2  \og{n)'^(T j y/n  and  apply  (soft)  thresholding  to 
the  wavelet  coefficients. 

3.  Invert  the  cascade  algorithm  to  get  the  denoised  signal  f[n]. 

Donoho  and  Johnstone’s  algorithm  offers  the  advantages  of  smoothness  and  adap¬ 
tation.  Wavelet  shrinkage  is  smooth  in  the  sense  that  the  denoised  estimate  /  has 
a  very  high  probability  of  being  as  smooth  as  the  original  signal  /,  in  a  variety  of 
smoothness  spaces  (Sobolev,  Holder,  etc.).  Wavelet  shrinkage  also  achieves  near- 
minimax  mean-square-error  among  possible  denoisings  of  /,  measured  over  a  wide 
range  of  smoothness  classes.  In  these  numerical  senses,  wavelet  shrinkage  is  superior 
to  other  smoothing  and  denoising  algorithms.  Heuristically,  wavelet  shrinkage  has 
the  advantage  of  not  adding  ‘ffiumps”  or  false  oscillations  in  the  process  of  removing 
noise,  because  of  the  local  and  smoothness-preserving  nature  of  the  wavelet  trans¬ 
form.  Wavelet  shrinkage  has  been  successfully  applied  to  SAR  imagery  as  a  method 
for  clutter  removal  [74].  It  is  natural  to  attempt  to  use  multiwavelets  as  the  transform 
for  a  wavelet  shrinkage  approach  to  denoising,  and  compare  the  results  with  scalar 
wavelet  shrinkage. 

We  implemented  Donoho’s  wavelet  shrinkage  algorithm  using  several  additional 
remarks  from  [74].  We  compared  the  performance  of  the  D4  scalar  wavelet  transform 
with  oversampled  and  critically  sampled  multi  wavelet  schemes.  In  the  oversampled 
scheme,  the  first  row  is  multiplied  by  a/2?  to  better  match  the  first  eigenvector  of 
the  GHM  system.  The  critically  sampled  scheme  obtains  two  input  rows  '^2,n 
from  a  single  row  of  data.  After  reconstruction  the  two  output  rows  Vi^ni  '^2,n  ^'I'e 
deapproximated  to  yield  the  output  signal  f[n].  Boundaries  are  handled  by  symmet¬ 
ric  data  extension  for  the  critically  sampled  (approximation/deapproximation)  and 
oversampled  schemes,  and  by  circular  periodization  for  D4. 
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Combining  these  various  techniques,  we  were  able  to  apply  multiwavelet  denoising 
to  imagery.  We  added  white  Gaussian  noise  to  the  Lenna  image,  and  applied  three 
wavelet  transforms  for  denoising  by  wavelet  shrinkage:  the  GHM  multi  wavelet  with 
approximation,  GHM  with  repeated  row,  and  the  Daubechies  4-tap  scalar  wavelet. 
The  experimental  results  are  shown  in  Table  4  and  the  resulting  images  in  Figure  47. 
The  GHM  multiwavelet  with  approximation  was  superior  to  D4  both  numerically  and 
subjectively;  the  approximation-based  preprocessing  seemed  to  reduce  the  Cartesian 
artifacts  present  in  the  scalar  wavelet  shrinkage.  This  can  be  seen,  for  example,  in 
the  facial  features  (eyes,  nose)  of  the  Lenna  images  shown  in  Figure  47.  The  GHM- 
repeated  row  scheme  suffered  because  we  had  to  repeat  rows  in  first  the  x  dimension 
and  then  in  the  y  dimension,  altering  the  correlations  of  the  data.  This  produces  the 
broad  stripes  in  the  image  denoised  with  the  repeated  row  scheme. 


Figure  47:  See  attached  page 


Noise 

GHM  with 
approximation 

GHM  with 
repeated  row 

£>4 

V-  error 

19.93 

6.87 

9.70 

8.09 

error 

24.98 

9.75 

12.6 

11.4 

Table  4:  Denoising  of  Lenna  image  via  wavelet-shrinkage. 


6.4.3  Transform-based  image  coding 

One  of  the  most  successful  applications  of  the  wavelet  transform  is  image  compression. 
A  transform-based  coder  operates  by  transforming  the  data  to  remove  redundancy, 
then  quantizing  the  transform  coefficients  (a  lossy  step),  and  finally  entropy  coding 
the  quantizer  output.  Because  of  their  energy  compaction  properties  and  correspon¬ 
dence  with  the  human  visual  system,  wavelet  representations  have  produced  superior 
objective  and  subjective  results  in  image  compression  [64],  [130],  [11],  [16].  Since  a 
wavelet  basis  consists  of  functions  with  short  support  (for  high  frequencies)  and  long 
support  (for  low  frequencies),  large  smooth  areas  of  an  image  may  be  represented 
with  very  few  bits,  and  detail  added  where  it  is  needed.  Multiwavelet  decompositions 
offer  all  of  these  traditional  advantages  of  wavelets,  as  well  as  the  combination  of  or¬ 
thogonality,  short  support,  and  symmetry.  The  short  support  of  multi  wavelet  filters 
limits  ringing  artifacts  due  to  subsequent  quantization.  Symmetry  of  the  filter  bank 
both  leads  to  efficient  boundary  handling  and  preserves  centers  of  mass,  lessening  the 
blurring  of  fine-scale  features.  Orthogonality  is  useful  because  it  means  that  rate- 
distortion  optimal  quantization  strategies  may  be  employed  in  the  transform  domain 
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Lenna  image  with  Gaussian  noise 
MSE  24.98 


Daubechies  4  scalar 
wavelet  denoising,  MSE  11.4 


Figure  47: 


GHM-with-approximation 
multiwavelet  denoising,  MSE  9.75 


GHM-repeated-row 
multiwavelet  denoising,  MSE  12.6 


comparison. 


and  still  lead  to  optimal  time-domain  quantization  (at  least  when  error  is  measured 
in  a  mean-square  sense).  Thus  it  is  natural  to  consider  the  use  of  multiwavelets  in  a 
transform-based  image  coder. 

We  compared  the  new  two-dimensional  multiwavelet  algorithms  with  a  scalar 
wavelet  in  a  production  image  coding  system.  Five  types  of  wavelet  transform  were 
used; 

•  £>4  scalar  wavelet 

•  Approximation/deapproximation  preprocessing  with  GHM  multiwavelets 

•  Adjacent  rows  input  with  GHM  multiwavelets 

•  Adjacent  rows  input  with  symmetric  pair 

•  Adjacent  rows  input  with  two  different  constrained  pairs 

Each  of  these  wavelet  transforms  was  followed  by  entropy-constrained  scalar  quanti¬ 
zation  and  entropy  coding.  We  made  the  assumption  that  the  histograms  of  subband 
(or  wavelet  transform  subblock)  coefficient  values  obeyed  a  Laplacian  distribution 
[64],  and  designed  a  uniform  scalar  quantizer.  The  quantizer  optimized  the  bit  allo¬ 
cation  among  the  different  subbands  by  using  an  operational  rate-distortion  approach 
(minimizing  the  functional  D  -f  \R)  [93].  We  then  entropy-coded  the  resulting  coef¬ 
ficient  streams  using  a  combination  of  zero-run-length  coding  and  adaptive  Huffman 
coding,  as  in  the  FBPs  Wavelet  Scalar  Quantization  standard  [7]. 


Compression  Ratio 

8:1 

16:1 

32:1 

64:1 

pSNR 

pSNR 

pSNR 

pSNR 

Daubechies  4 

35.6 

32.3 

29.3 

26.8 

GHM  with  appr./deappr. 

35.3 

31.8 

29.4 

27.1 

Adjacent  Row  Processing: 

GHM 

24.1 

21.3 

19.7 

18.4 

symmetric  pair 

31.1 

27.3 

24.0 

21.8 

constrained  pair 

32.4 

28.5 

25.1 

23.0 

constrained  pair  ^2 

31.9 

28.2 

25.0 

22.8 

Table  5:  Peak  SNRs  for  compression  of  Lenna. 

We  applied  these  different  wavelet  image  coders  to  the  Lenna  (NITF6)  image,  as 
well  as  a  geometric  test  pattern,  at  a  variety  of  compression  ratios.  The  results  are 
shown  in  Tables  5  and  6,  and  in  Figures  48  and  49.  On  Lenna,  the  GHM  multiwavelet 
with  approximation  mildly  outperformed  the  D4  scalar  wavelet  at  compression  ratios 
of  32:1  and  64:1.  The  images  in  Figure  48  show  that  the  GHM-approximation  scheme 
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preserves  more  texture  in  the  hat  and,  as  in  the  denoising  application,  produce  fewer 
Cartesian  artifacts  than  the  scalar  wavelet  scheme.  The  repeated-row  method  did 
not  work  well  on  Lenna.  However,  the  repeated-row  method  produced  the  best  com¬ 
pressions  of  the  test  pattern  image  (49)  at  intermediate  compression  ratios  (16:1 
and  32:1),  with  the  constrained  pair  #1  “CP-l”  outperforming  both  and  GHM 
with  approximation.  When  using  the  repeated  row  algorithm,  the  constrained  pairs 
significantly  outperformed  the  GHM  symmetric  multi  wavelet,  demonstrating  the  im¬ 
portance  of  the  eigenvector  constraints  used  in  definition  of  the  constrained  pairs.  A 
close  look  at  the  details  of  the  compressed/decompressed  test  patterns  shows  that 
the  CP-1  compression  did  a  better  job  of  preserving  the  checkerboard  and  “rang” 
over  a  shorter  distance  than  the  D4  compression.  The  alteration  of  the  checkerboard 
pattern  in  the  D4  compression  may  be  due  to  the  lack  of  linear-phase  symmetry  in 
the  wavelet  filters. 

These  preliminary  results  suggest  that  multiwavelets  are  worthy  of  further  inves¬ 
tigation  as  a  technique  for  image  compression.  Issues  to  address  include  the  design 
of  multiwavelets  with  symmetry  and  higher  order  of  approximation  than  the  GHM 
system,  the  role  of  eigenvector  constraints,  and  also  further  exploration  of  regular¬ 
ity  for  multiwavelets  [119].  One  might  also  apply  zerotree-coding  methods  [92]  in  a 
multiwavelet  context. 


Compression  Ratio 

8:1 

16:1 

32:1 

64:1 

pSNR 

pSNR 

pSNR 

pSNR 

Daubechies  4 

48.5 

31.4 

23.0 

19.8 

GHM  with  appr./deappr. 

52.4 

34.0 

18.3 

16.8 

Adjacent  Row  Processing: 

GHM 

29.8 

25.4 

20.1 

15.8 

symmetric  pair 

33.3 

28.3 

20.9 

16.5 

constrained  pair  1 

42.2 

32.3 

23.9 

19.0 

constrained  pair  2 

33.3 

30.2 

21.6 

17.6 

Table  6:  Peak  SNRs  for  compression  of  geometric  test  pattern. 


Figure  48:  See  attached  page 
Figure  49:  See  attached  page 
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Daubechies  4  Constrained-Pair  #1 

scalar  wavelet  compression  multiwavelet  compression 

64:1 ,  pSNR  26.8  64:1 ,  pSNR  23.0 


Figure  48:  Lenna  compression  comparison. 


Original  geometric  pattern 


Detail  of  original  pattern 
(corner  of  lower  left  checkerboard) 


Detail  of  CP-1  multiwavelet  Detail  of  D4  scalar  wavelet 

compression  (32:1,  pSNR  23.9)  compression  (32:1,  pSNR  23.0) 


Figure  49:  Pattern  compression  comparison. 


In  summary,  we  have  applied  the  new  mathematical  constructions  of  multiwavelets. 
Multiwavelets  offer  the  advantages  of  combining  symmetry,  orthogonality,  and  short 
support,  properties  not  mutually  achievable  with  scalar  2-band  wavelet  systems.  How¬ 
ever,  multiwavelets  differ  from  scalar  wavelet  systems  in  requiring  two  or  more  input 
streams  to  the  multiwavelet  filter  bank.  We  described  two  methods  (repeated  row 
and  approximation/deapproximation)  for  obtaining  such  a  vector  input  stream  from 
a  one-dimensional  signal.  We  developed  the  theory  of  symmetric  extension  for  multi¬ 
wavelet  filter  banks,  which  matches  nicely  with  approximation-based  preprocessing. 
We  then  applied  this  arsenal  of  techniques  to  two  basic  signal  processing  problems, 
denoising  via  thresholding  (wavelet  shrinkage),  and  data  compression.  After  devel¬ 
oping  the  approach  via  model  problems  in  one  dimension,  we  applied  the  various 
new  multiwavelet  approaches  to  the  processing  of  images,  frequently  obtaining  per¬ 
formance  superior  to  the  comparable  scalar  wavelet  transform.  These  results  suggest 
that  further  work  in  the  design  and  application  of  multiwavelets  to  signal  and  image 
processing  is  well  warranted. 


7  Technical  Conclusions 

In  conclusion,  this  project  explored  and  delineated  the  theory  of  a  broad  new  class  of 
mathematical  transforms,  the  rank  M  wavelets.  We  developed  a  production-quality 
software  package  (WaveTool)  for  the  design  and  implementation  of  these  transforms, 
and  applied  them  to  concrete  problems  in  signal  processing  and  communications.  In 
the  theoretical  realm,  we  discovered  complete  parametrizations  of  the  family  of  rank 
M  wavelets,  and  developed  and  applied  mathematical  tools  for  measuring  the  Sobolev 
smoothness  of  rank  M  wavelet  systems,  leading  to  surprising  asymptotic  regularity 
results.  We  also  developed  numerous  techniques  for  rank  M  wavelet  filter  design, 
including  those  based  on  approximation  (vanishing  moments)  and  on  smoothness  of 
the  iterated  filter  (regularity).  Methods  for  the  construction  of  full  rank  M  wavelet 
matrices  led  directly  to  fast  algorithms  for  computation,  particularly  in  the  case  of 
the  cosine-modulated  filter  banks. 

This  project  also  saw  the  design  and  completion  of  a  broad  and  flexible  software 
system,  “WaveTool”,  for  prototyping  wavelet  algorithms.  After  installation  at  a  num¬ 
ber  of  government,  academic,  and  industrial  beta  sites,  this  software  was  successfully 
turned  into  a  commercial  product.  Finally,  new  applications  of  rank  M  wavelets  and 
design  techniques  were  explored,  particularly  to  multicarrier  modulation  for  broad¬ 
band  communications.  This  latter  application  has  led  directly  to  a  chipset  product 
for  high-bitrate  communications  across  twisted-pair  copper  wire  and  hybrid  fiber-coax 
networks.  We  also  improved  our  leading-edge  wavelet  image  compression  algorithms, 
and  applied  wavelet-based  compression  to  sonar,  seismic,  and  multispectral  image 
data.  Finally,  we  produced  the  first  significant  applications  of  the  new  multiwavelet 


81 


techniques  to  signal  and  image  processing,  yielding  promising  results  in  denoising  and 
compression. 

Directions  for  further  research  that  were  suggested  by  this  work  include: 

•  The  systematic  application  of  wavelet-based  smoothness  measures  in  image 
quantization  and  compression.  Initial  steps  have  been  taken  in  [19]. 

•  Deeper  exploration  of  cosine-modulated  filter  bank  and  wavelet  structures,  in¬ 
cluding  filter  design  and  fast  algorithms  for  computation  [71],  [42]. 

•  Thorough  development  of  wavelet-based  systems  for  multicarrier  modulation, 
addressing  issues  such  as  equalization  and  system  latency  [47]. 

•  Deeper  exploration  of  multiwavelet  techniques  for  signal  processing,  incuding 
advances  in  multiwavelet  filter  design  and  preprocessing  algorithms  [129]. 

Of  course  this  is  only  a  partial  list;  many  other  topics  can  be  named. 

8  Participants 

Employees  at  Aware  who  participated  in  work  on  this  contract  included.  Dr.  Peter 
Niels  Heller,  Dr.  Howard  L.  Resnikoff,  Dr.  Michael  Tzannes,  Dr.  John  Weiss, 
Dr.  Edmund  Reiter,  Hemant  Singh,  Lev  Weisfeiler,  Vasily  Strela,  W.  Knox  Carey, 
Dr.  Richard  Gross,  Dr.  Stephen  DelMarco,  Karl  Jagler,  L.  Scott  Hills,  and  Anna 
Rounbehler. 

A  number  of  consultants,  each  of  them  experts  in  their  respective  fields,  were  em¬ 
ployed  as  well.  They  included:  Professor  R.  0.  Wells,  Jr.  of  Rice  University  (wavelet 
mathematics  and  computation).  Professor  P.  P.  Vaidyanathan  of  Cal  Tech  (multi¬ 
rate  signal  processing).  Professor  Truong  Q.  Nguyen  of  the  University  of  Wisconsin 
(design  of  wavelet  and  filter  bank  structures).  Dr.  Ramesh  Gopinath  of  Rice  Univer¬ 
sity  (wavelets  and  signal  processing),  and  Dr.  Sundar  Narasimhan  of  MIT  (graphics 
software  for  WaveTool). 


9  Publications 

Publications  arising  from  the  work  completed  under  this  contract  include  five  refereed 
journal  articles,  14  conference  papers,  and  two  book  chapters.  Specifically,  they 
include: 
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Journal  Articles: 


1.  P.  Steffen,  P.  N.  Heller,  R.  A.  Gopinath,  C.  S.  Burrus,  “Theory  of  Regular 
M-band  Wavelets,”  IEEE  Trans,  on  SP,  41  (1993),  pp.  3497-3511. 

2.  P.  N.  Heller,  “Rank  M  Wavelets  With  N  Vanishing  Moments,”  SIAM  J.  Matrix 
Analysis,  16  (1995),  pp.  502-519. 

3.  H.  L.  Resnikoff,  “Analytic  Representation  of  Compactly  Supported  Wavelets,” 
in  Festschrift  for  Hans  Bremermann,  BioSystems,  34  (1995),  pp.  259-272. 

4.  V.  Strela,  P.  N.  Heller,  G.  Strang,  P.  Topiwala,  and  C.  Heil,  “The  Application 
of  Multiwavelet  Filter  Banks  to  Image  Processing,”  to  appear,  IEEE  Trans,  on 
Image  Processing. 

5.  P.  N.  Heller  and  R.  0.  Wells,  Jr.,  “Sobolev  Regularity  for  Rank  M  Wavelets,” 
submitted  to  SIAM  J.  Math.  Analysis,  1996. 

Conference  Papers: 

1.  P.  N.  Heller  and  H.  L.  Resnikoff,  “Regular  M-band  wavelets  and  applications,” 
Proc.  IEEE  ICASSP,  Minneapolis,  1993. 

2.  C.  Bosman  and  E.  C.  Reiter,  “Seismic  data  compression  using  wavelet  trans¬ 
forms,”  SEG  Annual  Meeting  Extended  Abstracts,  pp.  1261-1264,  1993. 

3.  P.  N.  Heller  and  K.  Jagler,  “Wavelet  compression  of  multispectral  imagery,” 
in  Proc.  Industry  Workshop  —  Data  Compression  Conference,  Snowbird,  Utah, 
1994. 

4.  H.  L.  Resnikoff,  “Perfect  reconstruction  and  wavelet  matrix  windows  for  har¬ 
monic  analysis,”  in  Proc.  SPIE,  San  Diego,  CA,  1994. 

5.  P.  N.  Heller,  “Lagrange  M-th  band  filters  and  the  construction  of  smooth  M- 
band  wavelets,”  in  Proc.  lEEE-SP  Inti.  Symp.  on  Time-Frequency  and  Time- 
Scale  Analysis,  Philadelphia,  PA,  1994,  pp.  108-111. 

6.  E.  C.  Reiter  and  P.  N.  Heller,  “Wavelet  transform  based  compression  of  NMO- 
corrected  CDP  gathers,”  Society  of  Exploration  Geophysicists  64th  Annual 
Mtg.,  Los  Angeles,  1994.,  pp.  731-734. 

7.  M.  A.  Tzannes,  M.  C.  Tzannes,  J.  Proakis,  P.  N.  Heller,  “DMT  Systems, 
DWMT  systems,  and  digital  filter  banks,”  in  Proc.  IEEE  ICC,  New  Orleans, 
LA,  1994. 
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8.  P.  N.  Heller,  J.  M.  Shapiro,  and  R.  0.  Wells,  Jr.,  “Optimally  smooth  symmet¬ 
ric  quadrature  mirror  filters  for  image  coding,”  in  Proc.  SPIE  2491,  Wavelet 
applications  for  dual  use,  Orlando,  FL,  April,  1995. 

9.  P.  N.  Heller,  V.  Strela,  G.  Strang,  P.  Topiwala,  C.  Heil,  L.  S.  Hills,  “Multi¬ 
wavelet  filter  banks  for  data  compression,”  in  Proc.  IEEE  ISCAS  ’95,  Seattle, 
Washington. 

10.  P.  N.  Heller,  T.  Q.  Nguyen,  H.  Singh,  W.  K.  Carey,  “Linear-Phase  M-band 
wavelets  with  application  to  image  coding,”  in  Proc.  IEEE  ICASSP,  Detroit, 
MI,  1995. 

11.  S.  DelMarco,  P.  N.  Heller,  and  J.  Weiss,  “An  M-band,  2-dimensional  translation- 
invariant  wavelet  transform  and  applications,”  in  Proc.  IEEE  ICASSP,  Detroit, 
MI,  1995. 

12.  H.  Singh  and  P.  N.  Heller,  “WaveTool:  an  integrated  software  for  wavelet  and 
multirate  signal  processing,”  IEEE  Inti.  Conf.  on  Image  Proc.,  Washington, 
D.C.,  1995. 

13.  V.  Strela,  P.  N.  Heller,  G. Strang,  P.Topiwala  and  C.Heil,  “Application  of  mul¬ 
tiwavelets  to  signal  and  image  processing,”  in  UK  Symposium  on  Applications 
of  Time-Frequency  and  Time-Scale  Methods,  Warwick,  UK,  1995. 

14.  M.  Lang  and  P.  N.  Heller,  “The  Design  of  maximally  smooth  wavelets,”  IEEE 
ICASSP,  Atlanta,  GA,  1996. 

Book  Chapters: 

1.  P.  N.  Heller  and  R.  0.  Wells,  Jr.,  “The  Spectral  Theory  of  Multiresolution  Op¬ 
erators  and  Applications,”  in  Wavelets:  Theory,  Algorithms,  and  Applications, 
C.  K.  Chui,  L.  Montefusco,  L.  Puccio,  eds.,  Academic  Press,  San  Diego,  1994, 
pp.  13-32. 

2.  P.  N.  Heller,  “Tutorial  2.5:  Subband  and  Wavelet  Transforms;  Theory,  De¬ 
sign,  and  Applications:  Educational  Software,”  in  Microsystems  Technology  for 
Multimedia  Applications,  IEEE  ISCAS  Tutorial  volume,  1995. 

In  addition  a  number  of  invited  talks  were  given,  including: 

•  Nordic  Postgraduate  course  on  Wavelets  and  Filter  Banks,  Helsinki,  Finland, 
1994  -  P.  N.  Heller 

•  Sherman  Memorial  Lecture  at  Indiana  University,  Bloomington,  IN,  1994  -  H. 
L.  Resnikoff 
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•  Presentation  at  the  1993  Taormina,  Italy  conference  on  Wavelets  and  Applica¬ 
tions  -  R.  0.  Wells,  Jr. 

•  Hour  talk  at  Argonne  Workshop  on  Wavelets  and  Large-Scale  Image  Processing, 
Argonne,  IL,  1994  -  P.  N.  Heller 

10  Transition  of  Technology  to  Government  and 
Commercial  Uses 

10.1  WaveTool 

The  WaveTool  wavelet  and  multirate  design  and  algorithm  prototyping  software  was 
installed  at  the  following  government,  industrial,  and  academic  “beta  sites”: 

•  National  Security  Agency  (R5):  Adolf  Cusmariu  and  Mark  Marson 

•  National  Agency:  Laszlo  Fulop 

•  Analog  Devices,  Norwood,  MA 

•  General  Instrument,  Hatboro,  PA 

•  Massachusetts  Institute  of  Technology  (Math  Dept,  and  Civil  Engineering 
Dept.) 

•  California  Institute  of  Technology  (Electrical  Engineering  Dept.) 

•  Rice  University  (Electrical  and  Computer  Engineering  Dept,  and  Math  Dept.) 

•  University  of  California  at  Davis  (Electrical  Engineering  Dept.) 

•  Helsinki  University  of  Technology  (Electrical  Engineering  Dept.) 

Feedback  from  these  users  proved  very  useful  for  later  improvements  to  WaveTool, 
preceding  its  release  as  a  commercial  product  in  April  1995.  Since  then,  a  number  of 
copies  of  this  specialized  signal  processing  software  have  been  delivered  to  industry 
and  academic  customers  around  the  world. 

10.2  Other  commercial  products 

Work  done  on  this  contract  has  contributed  to  numerous  other  Aware  products  in 
both  last-mile  telecommunications  and  in  image  compression.  In  telecommunications, 
the  wavelet  filter  designs  and  computational  algorithms  of  section  4.4  as  applied  to 
multicarrier  modulation  (section  6.1)  have  been  designed  into  an  application-specific 
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integrated  circuit  (ASIC)  with  our  partner.  Analog  Devices.  This  chip,  the  AD6434, 
implements  a  384-tone  DWMT  algorithm  using  a  cosine-modulated  wavelet  filter  bank 
with  rank  M  =  384  and  genus  g  =  G.  It  is  intended  for  use  in  both  hybrid  fiber-coax 
and  twisted-pair  systems.  As  part  of  a  joint  project  with  DSC  Communications  Corp., 
the  AD6434  has  been  incorporated  into  DSC’s  MediaSpan  architecture  to  provide 
the  cable  telephony  subsystem.  Just  recently,  the  first  successful  phone  call  over  this 
DWMT-based  telephony  system  was  made.  Similar  DWMT  algorithms  are  being 
designed  into  an  Aware/ ADI  chipset  for  a  VDSL  system  (high-bitrate  communications 
over  twisted-pair  copper  lines  of  length  1000  to  3000  feet)  that  will  be  released  during 
1997. 

This  project’s  advances  in  wavelet  filter  design  and  image  compression  algorithms 
have  contributed  to  Aware’s  numerous  successful  products  for  wavelet-based  still  im¬ 
age  compression.  The  extension  of  wavelet  methods  to  seismic  data  compression  has 
led  directly  to  a  software  product,  SeisPact,  that  is  seeing  wide  use.  SeisPact  has 
been  widely  deployed  for  ship-to-shore  communication  of  seismic  exploration  results 
via  Inmarsat  (satellite)  links,  as  discussed  in  section  6.2.4,  SeisPact  is  also  being 
used  in  on-shore  seismic  data  processing  applications  to  minimize  data  storage  re¬ 
quirements.  Aware’s  AccuPress  software  for  compression  of  8-bit  and  24-bit  images 
derived  a  number  of  improvements  from  the  project  being  reported,  and  is  being  used 
for  a  variety  of  applications  from  medical  image  compression  to  multimedia.  In  addi¬ 
tion,  Aware  has  become  the  leading  vendor  of  compliance-certified  implementations 
of  the  FBI’s  Wavelet  Scalar  Quantization  [7]  wavelet-based  fingerprint  compression 
algorithm. 

10.3  Other  government  work 

Some  of  the  ideas  explored  during  this  contract  have  led  to  further  government  work 
as  well.  Aware  has  successfully  investigated  compression  of  one-dimensioal  acoustic 
data  via  a  Phase  I  SBIR  with  the  Navy  [26],  and  is  currently  working  on  Phase  II 
of  the  same  project,  intended  to  lead  to  a  compression  product  for  one-dimensional 
data. 

We  also  successfully  pursued  and  completed  a  Phase  I  SBIR  for  the  National  In¬ 
stitute  of  Standards  and  Technology,  the  “Wavelet  Image  Compression  Workbench,” 
[40]  that  took  many  of  the  ideas  from  this  contract  involving  choice  of  wavelet  basis 
and  multiresolution  tree  and  incorporated  them  into  a  software  plug-in  module  for 
Adobe  Photoshop.  This  led  directly  to  a  software  workbench  for  comparing  and  pro¬ 
totyping  wavelet  image  compression  algorithms,  and  to  a  compression  product  that 
we  are  selling  today. 

Finally,  new  applications  of  wavelet  techniques  to  communications  were  the  sub¬ 
ject  of  a  follow-on  government  agency  contract  with  Howard  L.  Resnikoff  [85]. 
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